{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weak Formulation\n",
    "\n",
    "Equation (b) is a vector equation, so really it is two or three equations involving multiple components of $\\vec{j}$. We want to work with a single scalar equation, allow for anisotropic physical properties, and potentially work with non-axis-aligned meshes — how do we do this?! We can use the **weak formulation** where we take the inner product ($\\int \\vec{a} \\cdot \\vec{b} dv$) of the equation with a generic face function, $\\vec{f}$. This reduces requirements of differentiability on the original equation and also allows us to consider tensor anisotropy or curvilinear meshes (Haber 2014).\n",
    "\n",
    "## Scalar equations only, please\n",
    "\n",
    "In Figure 5, we visually walk through the discretization of equation (b). On the left hand side, a dot product requires a *single* cartesian vector, $[\\mathbf{j_x, j_y}]$. However, we have a $j$ defined on each face (2 $j_x$ and 2 $j_y$ in 2D!). There are many different ways to evaluate this inner product: we could approximate the integral using trapezoidal, midpoint or higher order approximations. A simple method is to break the integral into four sections (or 8 in 3D) and apply the midpoint rule for each section using the closest $\\mathbf{j}$ components to compose a cartesian vector. A $\\mathbf{P}_i$ matrix (size 2 × 4) is used to pick out the appropriate faces and compose the corresponding vector (these matrices are shown with colors corresponding to the appropriate face in the figure). On the right hand side, we use a vector identity to integrate by parts. The second term will cancel over the entire mesh (as the normals of adjacent cell faces point in opposite directions) and $\\phi$ on mesh boundary faces are zero by the Dirichlet boundary condition. This leaves us with the  divergence, [which we already know how to do](divergence.ipynb)!\n",
    "\n",
    "<img src=\"images/WeakFormulation.png\" style=\"width:75%;text-align:center\">\n",
    "\n",
    "<p style=\"text-align:center;\">\n",
    "**Figure 5**: Discretization using the weak formulation and inner products.\n",
    "</p>\n",
    "\n",
    "The final step is to recognize that, now discretized, we can cancel the general face function $\\mathbf{f}$ and transpose the result (for convention's sake):\n",
    "\n",
    "$$\n",
    "\\frac{1}{4}\\sum_{i=1}^{4} \\mathbf{P}_i^\\top \\sqrt{v} \\boldsymbol{\\Sigma}^{-1} \\sqrt{v} \\mathbf{P}_i \\mathbf{j} = \\mathbf{D}^\\top v \\phi\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation\n",
    "\n",
    "We will start by using a SimPEG Mesh that only has a single cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from SimPEG import Mesh\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEw9JREFUeJzt3W2MXOV5xvHrchFKeCkuaeNUNNiW80KhSbbEMVhV5eUD\nBaK2bhvCS01h06hBieBDpEjQKsQ4pFIoKKmAIEKEkhg5Ii6GxCRAcYsPyNRsTIgTQ6BxqNcBAm4l\nWAkH1Fr47ocZm81yZl/OPnueM2f+P2nkOTPHu3eePN6Lue85s44IAQAwWwtyFwAA6E8ECACgEgIE\nAFAJAQIAqIQAAQBUQoAAACrJHiC2b7O9z/ZPejy/yva47ce7t8/WXSMA4M2OyF2ApK9LulHS+inO\neTgi/rymegAAM5D9FUhEbJP08jSnuY5aAAAzlz1AZmil7Z22v2/75NzFAACa0cKazg8lnRgRr9o+\nR9J3JL0nc00AMPAaHyARsX/C/fts32z7+Ih4afK5tvlgLwCYpYioNCZoSgvL6jHnsL1owv0VklwW\nHodEBLcEt7Vr12avoU031pP1bOptLrK/ArH9LUnDkt5m+xeS1ko6UlJExK2SzrX9SUkHJL0m6fxc\ntQ6SsbGx3CW0CuuZFuvZDNkDJCL+eprnvyLpKzWVAwCYoaa0sNAwIyMjuUtoFdYzLdazGTzXHliT\n2I42/e8BgPlmW9HnQ3Q0TFEUuUtoFdYzLdazGQgQAEAltLAAYIDRwgIA1I4AQSl6zGmxnmmxns1A\ngAAAKmEGAgADjBkIAKB2BAhK0WNOi/VMi/VsBgIEAFAJMxAAGGDMQAAAtSNAUIoec1qsZ1qsZzMQ\nIACASpiBAMAAYwYCAKgdAYJS9JjTYj3TYj2bgQABAFTCDAQABhgzEABA7QgQlKLHnBbrmRbr2QwE\nCACgEmYgADDAmIEAAGpHgKAUPea0WM+0WM9mIEAAAJUwAwGAAcYMBABQOwIEpegxp8V6psV6NgMB\nAgCohBkIAAwwZiAAgNoRIChFjzkt1jMt1rMZCBAAQCXMQABggDEDAQDUjgBBKXrMabGeabGezUCA\nAAAqYQYCAAOMGQgAoHYECErRY06L9UyL9WwGAgQAUAkzEAAYYMxAAAC1I0BQih5zWqxnWqxnM2QP\nENu32d5n+ydTnHOD7d22d9oeqrM+AEC57AEi6euSzur1pO1zJC2LiHdLulTSLXUVNsiGh4dzl9Aa\nEaH7H7pfzOfSYX82Q/YAiYhtkl6e4pTVktZ3zx2VdJztRXXUBqSw6Z5NuvnBm3XX9+7KXQqQVPYA\nmYETJD074fj57mOYR/SY04gIXX/79Xpl6Su6bv11vApJhP3ZDEfkLiC1kZERLVmyRJK0cOFCDQ0N\nHX65e2jTcTz98RlnSFLnWBru/slxteNHJRUalbRgwUMNqKcdx1u3do6b8O+ln44P3R8bG9NcNeI6\nENuLJd0TEe8vee4WSVsj4tvd46clrYqIfSXnch1IIrbEUs5NRGjleSs1esqoZEkhnfbkadq+cbvs\nSm+7Rxf7M502XAfi7q3MZkkXS5Lt0yWNl4UH0DSb7tmkXcfuemNnW9p1zC5mIWiN7C0s299S5zXp\n22z/QtJaSUdKioi4NSLutf1h2z+X9CtJH8tX7SAp9EbrAFU88tgjWv76cnmPNf7iuBa+Y6EiQtt2\nbNNH/uwjucvrc4XYn/k1ooWVCi2sdOxCEcO5y2iNoih462lC7M905tLCIkBQih4zmoz9mU4bZiAA\ngD5DgKCHIncBrcJ1C6kVuQuACBAAQEXMQFCKHjOajP2ZDjMQAEDtCBD0UOQuoFWYgaRW5C4AIkAA\nABUxA0EpesxoMvZnOsxAAAC1I0DQQ5G7gFZhBpJakbsAiAABAFTEDASl6DGjydif6TADAQDUjgBB\nD0XuAlqFGUhqRe4CIAIEAFARMxCUoseMJmN/psMMBABQOwIEPRS5C2gVZiCpFbkLgAgQAEBFzEBQ\nih4zmoz9mQ4zEABA7QgQ9FDkLqBVmIGkVuQuACJAAAAVMQNBKXrMaDL2ZzrMQAAAtSNA0EORu4BW\nYQaSWpG7AIgAAQBUxAwEpegxo8nYn+kwAwEA1I4AQQ9F7gJahRlIakXuAiACBABQETMQlKLHjCZj\nf6bDDAQAUDsCBD0UuQtoFWYgqRW5C4AIEABARcxAUIoeM5qM/ZkOMxAAQO0IEPRQ5C6gVZiBpFbk\nLgAiQAAAFTEDQSl6zGgy9mc6zEAAALUjQNBDkbuAVmEGklqRuwCIAAEAVMQMBKXoMaPJ2J/pMAMB\nANSOAEEPRe4CWoUZSGpF7gKgBgSI7bNtP237Z7avKHl+le1x2493b5/NUScA4NdNGSC2f9P2spLH\n35/im9teIOkmSWdJOkXShbZPKjn14Yg4tXv7QorvjTf79Oc+rTs336nOHGk4dzmtsHfPHq276CJt\nXbdO6y66SHv37MldUt9ifzZPzwCxfZ6kpyVtsv2k7Q9NePobib7/Ckm7I2JvRByQdIek1WXlJPp+\nmMLjex/XJXdfopUfXSlJ4g0Jc7N3zx7deOaZ+syGDVpXFPrMhg268cwzCZGK2J/NM9UrkH+Q9MGI\nGJL0MUm32/7L7nOpfqCfIOnZCcfPdR+bbKXtnba/b/vkRN8bk9jWq0te1egfjEoqtPKjKyf8Fx9m\n6xtXXaV1zzyjo9Xp2B8tad0zz+gbV12Vt7A+xf5sniOmeO43IuIFSYqIH9g+Q9L3bL9TUp3/j/1Q\n0okR8artcyR9R9J7ep08MjKiJUuWSJIWLlyooaEhDQ8PS3pjkMlx+fH4i+OSJC3t/DF6zKguvupi\nPfajx/TFtV/MXl+/Hf/XE09oh95othTdPw/+8peNqK/fjtmfaY4P3R8bG9OcRUTpTdJ/SFo26bFj\nJf27pP/t9fdmc5N0uqT7JxxfKemKaf7OHknH93guUN2qS1aFrlZorUKKOO3c0+LOzXfGwYMHc5fW\nl65esyb2dy5XOHzbL8XVa9bkLq0vsT/nR/fnZqWf4VO9AvmkpAW2T46In3Z/Or9i+2xJF8w9uiRJ\nOyS9y/ZiSS90v+6FE0+wvSgi9nXvr1Dn4seXEn1/TBAROmrsKL1v//s0Kmn7xu2yGT9VNXLNNVr7\n6KOH21i/krR22TJdfs01uUvrS+zP5uk5A4mIH0fEbkkbbV/hjrdK+pKkT6X45hHxuqTLJD0g6UlJ\nd0TEU7Yvtf2J7mnn2n7C9o8k/bOk81N8b7zZqYtP1fq/Wq/tG7dLKvjHOUeLly7V5Vu26Po1a3TJ\n0JCuX7NGl2/ZosVLl+YurS+xP5tn2o8ysX20pGslfVCdFtYGSddGxMH5L292+CiTdOxCEcO5y2iN\noigO96Ixd+zPdOb7o0wOSHpN0lslvUXSniaGB1Ibzl1AqxAeqQ3nLgCaWYDsUCdAPiTpj9W52O9f\n5rUqAEDjzSRAPh4Rn4uIAxHxQkSslrR5vgtDbkXuAlpl4lsokUKRuwBoBgESEY+VPHb7/JQDAOgX\n/D4QlOL3LaDJ2J/p8PtAAAC1I0DQQ5G7gFZhBpJakbsAiAABAFTEDASl6DGjydif6TADAQDUjgBB\nD0XuAlqFGUhqRe4CIAIEAFARMxCUoseMJmN/psMMBABQOwIEPRS5C2gVZiCpFbkLgAgQAEBFzEBQ\nih4zmoz9mQ4zEABA7QgQ9FDkLqBVmIGkVuQuACJAAAAVMQNBKXrMaDL2ZzrMQAAAtSNA0EORu4BW\nYQaSWpG7AIgAAQBUxAwEpegxo8nYn+kwAwEA1I4AQQ9F7gJahRlIakXuAiACBABQETMQlKLHjCZj\nf6bDDAQAUDsCBD0UuQtoFWYgqRW5C4AIEABARcxAUIoeM5qM/ZkOMxAAQO0IEPRQ5C6gVZiBpFbk\nLgAiQAAAFTEDQSl6zGgy9mc6zEAAALUjQNBDkbuAVmEGklqRuwCIAAEAVMQMBKXoMaPJ2J/pMAMB\nANSOAEEPRe4CWoUZSGpF7gIgAgQAUBEzEJSix4wmY3+mwwwEAFA7AgQ9FLkLaBVmIKkVuQuAGhAg\nts+2/bTtn9m+osc5N9jebXun7aG6awQAvFnWALG9QNJNks6SdIqkC22fNOmccyQti4h3S7pU0i21\nFzqQhnMX0BoRofsful/M51Iazl0AlP8VyApJuyNib0QckHSHpNWTzlktab0kRcSopONsL6q3TKC6\nTfds0s0P3qy7vndX7lKApHIHyAmSnp1w/Fz3sanOeb7kHCRX5C6gFSJC199+vV5Z+oquW38dr0KS\nKXIXAElH5C4gtZGRES1ZskSStHDhQg0NDWl4eFjSG4NMjmd2bHeO32gXcFzt+FFJhUYlLVjwUAPq\nacdx7n8f/Xp86P7Y2JjmKut1ILZPl3R1RJzdPb5SUkTEtRPOuUXS1oj4dvf4aUmrImJfydfjOhA0\nRkRo5XkrNXrKqGRJIZ325GnavnG77EpvuweS6+frQHZIepftxbaPlHSBpM2Tztks6WLpcOCMl4UH\n0DSb7tmkXcfu6oSHJFnadcwuZiFojawBEhGvS7pM0gOSnpR0R0Q8ZftS25/onnOvpD22fy7pq5I+\nla3gAcJ1C3P3yGOPaPnry7Vqzyp9YPsHtGrPKi0/uFzbdmzLXVrfY382Q/YZSETcL+m9kx776qTj\ny2otCkjgy5//8uH7RVEc7kUDbcFnYQHAAOvnGQgAoE8RIChFjzkt1jMt1rMZCBAAQCXMQABggDED\nAQDUjgBBKXrMabGeabGezUCAAAAqYQYCAAOMGQgAoHYECErRY06L9UyL9WwGAgQAUAkzEAAYYMxA\nAAC1I0BQih5zWqxnWqxnMxAgAIBKmIEAwABjBgIAqB0BglL0mNNiPdNiPZuBAAEAVMIMBAAGGDMQ\nAEDtCBCUosecFuuZFuvZDAQIAKASZiAAMMCYgQAAakeAoBQ95rRYz7RYz2YgQAAAlTADAYABxgwE\nAFA7AgSl6DGnxXqmxXo2AwECAKiEGQgADDBmIACA2hEgKEWPOS3WMy3WsxkIEABAJcxAAGCAMQMB\nANSOAEEpesxpsZ5psZ7NQIAAACphBgIAA4wZCACgdgQIStFjTov1TIv1bAYCBABQCTMQABhgzEAA\nALXLFiC2f8v2A7b/0/a/2j6ux3ljtn9s+0e2f1B3nYOKHnNarGdarGcz5HwFcqWkf4uI90p6UNLf\n9zjvoKThiPjDiFhRW3UDbufOnblLaBXWMy3WsxlyBshqSd/s3v+mpL/ocZ5Fq6124+PjuUtoFdYz\nLdazGXL+YH57ROyTpIh4UdLbe5wXkrbY3mH772qrDgAwpSPm84vb3iJp0cSH1AmEz5ac3uvtU38U\nES/Y/h11guSpiNiWuFRMMjY2lruEVmE902I9myHb23htP6XObGOf7XdI2hoRvz/N31kr6ZWI+FKP\n53kPLwDMUtW38c7rK5BpbJY0IulaSZdI+u7kE2wfJWlBROy3fbSkP5G0rtcXrLoIAIDZy/kK5HhJ\nGyW9U9JeSedFxLjt35X0tYj4U9tLJd2tTnvrCEkbIuKLWQoGAPyaVl2JDgCoT9++PZYLEdOwfbbt\np23/zPYVPc65wfZu2zttD9VdYz+Zbj1tr7I9bvvx7q3sDSWQZPs22/ts/2SKc9ibMzTdelbZm30b\nIOJCxDmzvUDSTZLOknSKpAttnzTpnHMkLYuId0u6VNIttRfaJ2aynl0PR8Sp3dsXai2yv3xdnbUs\nxd6ctSnXs2tWe7OfA4QLEeduhaTdEbE3Ig5IukOddZ1otaT1khQRo5KOs71IKDOT9ZQ6exLT6L5d\n/+UpTmFvzsIM1lOa5d7s5x+sXIg4dydIenbC8XPdx6Y65/mSc9Axk/WUpJXdlsv3bZ9cT2mtxN5M\nb1Z7M+fbeKfFhYhooR9KOjEiXu22YL4j6T2ZawKkCnuz0QESEWf2eq47DFo04ULE/+7xNV7o/vk/\ntu9Wp81AgHQ8L+nECce/131s8jnvnOYcdEy7nhGxf8L9+2zfbPv4iHipphrbhL2ZUJW92c8trEMX\nIkpTXIho+5ju/UMXIj5RV4F9YIekd9lebPtISReos64TbZZ0sSTZPl3S+KHWId5k2vWc2KO3vUKd\nt9ITHr1Zvfvy7M3Z67meVfZmo1+BTONaSRtt/626FyJK0sQLEdVpf93d/YiTQxciPpCr4KaJiNdt\nXybpAXX+Y+K2iHjK9qWdp+PWiLjX9odt/1zSryR9LGfNTTaT9ZR0ru1PSjog6TVJ5+eruNlsf0vS\nsKS32f6FpLWSjhR7s5Lp1lMV9iYXEgIAKunnFhYAICMCBABQCQECAKiEAAEAVEKAAAAqIUAAAJUQ\nIEANbN9n+2Xbky/UBPoWAQLU458kXZS7CCAlAgRIyPby7i8wO9L20bafsH1yRGyVtH/aLwD0kX7+\nKBOgcSLiMdvflfSPkt4q6faI+GnmsoB5QYAA6V2jzgcrvibp8sy1APOGFhaQ3m9LOkbSsZLekrkW\nYN4QIEB6t6jzS882qDM8P2SqjyYH+g4tLCAh238j6f8i4g7bCyQ9YntY0uclvVfSMd2P0v54RGzJ\nWCowZ3ycOwCgElpYAIBKCBAAQCUECACgEgIEAFAJAQIAqIQAAQBUQoAAACohQAAAlfw/RamZT0QY\nxqUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d6e2350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "mesh = Mesh.TensorMesh([1,1])\n",
    "mesh.plotGrid(ax=ax,centers=True,faces=True)\n",
    "ax.set_xlim(-0.5,1.5)\n",
    "ax.set_ylim(-0.5,1.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The figure above is of a single cell, which in 2D has four faces (the green triangles). Two of these are pointed in the x direction and two in the y direction.\n",
    "\n",
    "For this single cell, we need our projection matrices to pick out four cartesian vectors that we will then use in evaluating the full integral.\n",
    "\n",
    "Instead of creating the code from scratch in this notebook, we will leverage some of the internal functions that are in SimPEG. The code is straight forward, but can get a little bit intricate when working with multiple dimensions and indexing in a generalized way (i.e. works for all of the meshes).\n",
    "\n",
    "If you ever want to see the exact code that is executed in a function, you can use two question marks (one for the docs only). For example, execute **`mesh._getFacePxx??`** and Jupyter will show you more than you ever wanted. :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# mesh._getFacePxx??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function creates the P matrices that we need. The matrix just pick out the appropriate faces and are used later on to evaluate our full integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0.  0.]\n",
      " [ 0.  0.  1.  0.]]\n"
     ]
    }
   ],
   "source": [
    "# face-x minus side and the face-y minus side\n",
    "P1 = mesh._getFacePxx()('fXm','fYm')\n",
    "print( P1.todense() )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create and plot all of the matrices that were in Figure 5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA10AAACGCAYAAAAmauAuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC5ZJREFUeJzt3V+IbedZBvDnjaf2BBXBemFtaKvUsbaIGtuIVOioSNIW\n1AsFq9TQmyD4JyCIpTdnTvWiXongVTCtKNQKRagXUi02c0gttMWkNv3nKYJpFQ0IEdH2aKifF2ef\nMGdm9sx4sr+11rf37wcDM3uGvR6+vZ418+619p5qrQUAAIA+7po7AAAAwDYzdAEAAHRk6AIAAOjI\n0AUAANCRoQsAAKAjQxcAAEBHhi4AAICODF0dVdWDVXVl7hy7wFrvLo/9dKz1bvP4T8da7y6P/XSm\nXuvFDV1V9bqq+ruq+vqq+oaq+kxVvWbuXKe5YNbF/vfpqrpaVQ8f+fq3q+pX58y0zgWzLnatl0bP\npqNnu03XpqNru0vPpqNnLyBPa8t7XKvqXUnuXn18ubX2OzNHWuusrFX1YJJXttauzpXvLFX1iiR/\n1lr7waqqJF9M8rrW2r/PHO2ENVlf31p7dvX9Ra/1EunZNPQMXZuGru02PZuGnt25S1Ns5A78VpJP\nJvlqkkVOz0fclrWqviXJX+fm5PySJC+qqp9eff221tpnZ0t6TGvt6ar6t6r6viTfluSJJZYmOT1r\nkqqqJzPAWi+Unk1Az4iuTULXdp6eTUDP7txSz3S9NMnjSW4kuS/JO5O8JUlrrd07Z7bjjmV9fWvt\nq0e+92CSV7TW3jVXvvNU1c8meUNu7ox/mORHsty1vi1ra+1DR763+LVeGj2bjp7tNl2bjq7tLj2b\njp7dYZaFDl0fTPInSb4jybe31hb7jMVZWQcpzouSPJWbZz2/qy1xh1g5K+sIa700ejYdPdttujYd\nXdtdejYdPbszi7u8sKreluR/Wmvvr6q7kvxNVe231g5njnbCSFnXaa09V1WPJXl2yaVJxsq6dCPt\nuyNlXWekfXekrCMYaf8dKes6I+2/I2VdupH23ZGyrjPSvrukrIs808V0VoX/2yQ/01r7h7nznGWk\nrHDUSPvuSFnhuJH235GywlEj7btLyrq4t4xnOlX1Pbn5Ti4fnntHPM9IWeGokfbdkbLCcSPtvyNl\nhaNG2neXltWZLgAAgI6c6QIAAOjI0AUAANCRoQsAAKCjjb1lfFV5cRhbqbVWc2e4Rc/YVnoG/S2p\nZ4musb1O69pGz3S11jb+ceXKlS732/tjxNwjZu6de4lGXEeZx8+tZ8tfx1Eff5mnybxUo63jqI+/\nzNPlXsflhQAAAB0ZugAAADpa/NC1v78/d4Q7MmLuETMn4+ZemhHXccTMyZi5R8y8RKOu44i5Zd5d\no67jiLlHzJzMk3tj/xy5qtqm7guWoqrSFvTCYz1jG+kZ9Le0niW6xnZa17XFn+kCAAAYmaELAACg\no439n66zPPTQu3P9+o0pNnWbvb3LeeSRd0y+XeDiHB8YjX0WWMfxgXUmGbquX7+Ra9cOptjUMXNs\nE/j/cHxgNPZZYB3HB9ZxeSEAAEBHhi4AAICODF0AAAAdGboAAAA6MnQBAAB0ZOgCAADoyNAFAADQ\nkaELAACgowsNXVX1QFV9oaquV9Vv9g4Fu0jPYBq6Bv3pGdzu3KGrqu5K8vtJ7k/y2iRvrapX9w4G\nu0TPYBq6Bv3pGZx0kTNd9yX5Ymvt6dbac0nen+Sn+saCnaNnMA1dg/70DI65yND1siRfPvL1P61u\nAzZHz2Aaugb96Rkc4400AAAAOrp0gZ/55yQvP/L1PavbTjg4OHj+8/39/ezv77+AaDC9w8PDHB4e\nzrFpPWNnzNiz5IJd0zNGN0LPEl1jfBftWrXWzv6Bqq9L8vdJfjzJvyT5RJK3ttY+f+zn2rr72t8/\nyLVrBxfJvVFvfONBDg+n3y7bo6rSWqsJtvOCezYqxwem6tlqW+d27bye2WcZ0dJ6tvo5v9M2xPFh\nOdZ17dwzXa21r1XVryT5q9y8HPHR46UBXhg9g2noGvSnZ3DSRS4vTGvtQ0m+u3MW2Gl6BtPQNehP\nz+B23kgDAACgI0MXAABAR4YuAACAjgxdAAAAHRm6AAAAOjJ0AQAAdGToAgAA6MjQBQAA0JGhCwAA\noCNDFwAAQEeGLgAAgI4uTbGRvb3LSQ6m2NQp2wWWzPGB0dhngXUcH1inWmubuaOqtqn7gqWoqrTW\nau4ct+gZ20jPoL+l9SzRNbbTuq65vBAAAKAjQxcAAEBHhi4AAICODF0AAAAdGboAAAA6MnQBAAB0\nZOgCAADoyNAFAADQkaELAACgI0MXAABAR4YuAACAjgxdAAAAHRm6AAAAOro0dwA266GH3p3r129M\nvt29vct55JF3TL5dmIuuAadxbIBpjNY1Q9eWuX79Rq5dO5hhy3NsE+aja8BpHBtgGqN1zeWFAAAA\nHRm6AAAAOjJ0AQAAdGToAgAA6MjQBQAA0JGhCwAAoCNDFwAAQEeGLgAAgI7OHbqq6tGqeqaqPj1F\nINhVugb96Rn0p2dw0kXOdL03yf29gwC6BhPQM+hPz+CYc4eu1tpHkzw7QRbYaboG/ekZ9KdncJLX\ndAEAAHRk6AIAAOjo0ibv7ODg4PnP9/f3s7+/v8m7h+4ODw9zeHg4d4wz6Rmj0zPob4SeJbrG+C7a\ntYsOXbX6ONPR4sCIjh/wr169OnWEc7umZ4xOz6C/EXqW6Brju2jXLvKW8e9L8rEke1X1pap6+4Yy\nAkfoGvSnZ9CfnsFJ557paq39/BRBYNfpGvSnZ9CfnsFJ3kgDAACgI0MXAABAR4YuAACAjgxdAAAA\nHRm6AAAAOjJ0AQAAdGToAgAA6MjQBQAA0JGhCwAAoCNDFwAAQEeGLgAAgI4uzR2Azdrbu5zkYKbt\nwu7QNeA0jg0wjdG6Vq21jQSoqrap+4KlqKq01mruHLfoGdtIz6C/pfUs0TW207quubwQAACgI0MX\nAABAR4sfug4PD+eOcEdGzD1i5mTc3Esz4jqOmDkZM/eImZdo1HUcMbfMu2vUdRwx94iZk3lyG7o6\nGTH3iJmTcXMvzYjrOGLmZMzcI2ZeolHXccTcMu+uUddxxNwjZk4MXQAAAFvH0AUAANDRRt8yfiN3\nBAuzpLfY1TO2lZ5Bf0vqWaJrbK/TuraxoQsAAICTXF4IAADQkaELAACgI0MXAABAR4auBauqr1XV\nE1X1VFX9aVVdXt3+aFU9U1WfnjsjbIPTulZV91TVR6rqs6vbf23unDCyNT17cVV9vKqeXN1+Ze6c\nMLJ1fzuuvnfX6nt/PmfGXWXoWrb/aq3d21r73iTPJfml1e3vTXL/fLFg65zWteeS/Hpr7bVJfjjJ\nL1fVq+cMCYM70bPW2n8n+dHW2g8k+f4kb6qq+2ZNCWNb97djkjyc5HPzxMLQNY7Hk7wqSVprH03y\n7LxxYGs9nuRVrbVnWmufSpLW2n8m+XySl82aDLbH0d9pX1nd9uIkl5J4W2XYjOd7VlX3JHlzkj+Y\nNdEOM3QtWyVJVV1K8qYkT80bB7bWmV2rqlfm5rPwH586GGyRU3u2uuTpyST/muTDrbVPzhcRhne8\nZ7deivK7SX4jntSYjaFr2e6uqieSfCLJ00kenTkPbKu1Xauqb0zygSQPr854AXfm1J611v53dXnh\nPUl+qKpeM2NGGN3Rnv1jkvdU1VuS3Lp6o1YfTOzS3AE401daa/fOHQJ2wKldWz1T+IEkf9xa++D0\nsWCrnPk7rbX2H1X1WJIH4nUncKdO9Kyq3pDkJ6vqzUnuTvJNVfVHrbVfnCXhjnKma9nOeibCMxWw\nOeu69J4kn2ut/d6UYWBLnehZVX1rVX3z6vO7k/xEki9MHQy2yImetdbe2Vp7eWvtO5P8XJKPGLim\nZ+hatlOvu62q9yX5WJK9qvpSVb192liwdU50bfXM4C8k+bHV21k/UVUPTB8NtsZpv9NemuSxqvpU\nbr5m8i9ba38xbSzYKl6ztVDVmscGAACgF2e6AAAAOjJ0AQAAdGToAgAA6MjQBQAA0JGhCwAAoCND\nFwAAQEeGLgAAgI4MXQAAAB39HybJCU4iB2NNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c1bd110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P1 = mesh._getFacePxx()('fXm','fYm')\n",
    "P2 = mesh._getFacePxx()('fXp','fYm')\n",
    "P3 = mesh._getFacePxx()('fXm','fYp')\n",
    "P4 = mesh._getFacePxx()('fXp','fYp')\n",
    "fig, ax = plt.subplots(1,4, figsize=(15,3))\n",
    "def plot_projection(ii, ax, P):\n",
    "    ax.spy(P, ms=30)\n",
    "    if P.shape[1]==4:\n",
    "        ax.set_xticks(range(4))\n",
    "        ax.set_xticklabels(('x-','x+','y-','y+'))\n",
    "    ax.set_xlabel('P{}'.format(ii+1))\n",
    "map(plot_projection, range(4), ax, (P1, P2, P3, P4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This assumes that we unwrapped our face vectors in certain way and they are ordered [xfaces, yfaces]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# All the cells!\n",
    "\n",
    "When we move to more cells the visualization of the matrix is not as clear, so instead we will show a connection from each cell center to the x and y faces that are selected by the P matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mesh = Mesh.TensorMesh([3,4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABDEAAAGrCAYAAADZzwwPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+wZVd53vnnbcsEJNpoPFQxkWyrNU1sAhiEIqsll6y+\n7agK4bEjB4EL3A00nhTUzKA4FFQkZqIf7TilEDSVmOCJjYdIIwqXCrVUkWBsD3LSh4w0QhEVGrdl\niITULYOg7NgxDpKcikyv+eOe7vtDd22du895z3rfvb+fqlv07b6oHz17rbdPb521r5VSBAAAAAAA\nEN2O1gEAAAAAAABmwU0MAAAAAACQAjcxAAAAAABACtzEAAAAAAAAKXATAwAAAAAApMBNDAAAAAAA\nkAI3MQAAAAAAQArcxEBvZnbCzJ41s/9sZt8ys39pZmeZ2VvN7AEze8bM/k3rnADa6pgVHzGzR83s\nz83sD8zsHa2zAmirY1582Mz+cDovjpvZda2zAmhni1lxq5mdue7X/xsz+49m9m9b5oQPbmJgHkXS\n/1BK+T5JF0r6MUn/m6Q/lfRPJd3cMBuAOLaaFf9A0tOSfrqU8jJJByX9ipld0iwlgAhqry0+IenV\n03nx45IOmNnPtosJoLHNs+Iirb62OOXDkh5pEQz+uImBeZkklVK+Jem3Jb22lPJvSimHJX2raTIA\nkWyeFa8ppRwqpTw6/fl/J+n/lXRpu4gAgtjqtcWjpZSnp7++Q9JJSa9slA9ADM+bFZJkZj8u6TWS\nbm0XDZ64iYGFMLMflPRTkv596ywA4qrNCjN7iVb/iyv/1QSApOfPCzO71sy+I+nrks6U9JsN4wEI\nYv2sMLMdkv65pPe1TQVPZ7QOgPT+lZn9paQ/l/RZcYQEwNZeaFb8mqQvlVI+t/RkAKLZcl6UUj4s\n6cNm9npJPzv9dQDjtdWs+EVJD5ZSvmRmr2uaDm64iYF5XVVKOdI6BIDwqrPCzD4i6dWS9i03EoCg\nOl9blFK+bGZXSvolSR9YXiwAwWyYFWb2VyVdo9VnZEjT4yYYHm5iYF4MBwCz2HJWmNkhSW+UdPm6\n8+4Axm2W1xZnSPrvvYMACG3zrLhY0n8n6Q/MzCS9RNJLzOybks4tpZRlB4QPnomBhTOzHWb2VyR9\nr6TvMbO/YmbcMAOwgZl9SNLbJV1RSvl26zwAYrJV7zGzs6efXyzpf5H0u22TAQjmtyTtknSBpNdL\nukGrz9R5PTcwhoWbGJhHbRi8Q9JfSPpVSZdJelbSx5cVCkA4tVnxjyT9oKSvmdl3pt/r/bol5gIQ\nT21e/G2tzor/LOl2Sb9SSvnV5cUCEMzzZkUp5blSyh+f+tDqszKeK6X8x+XHgyfzvCllZp+Q9NOS\n/qiUsuWDVczso5LeJOkZSQdLKUfdAgEIi3kBYBbMCgCzYFYAw+X9ToxbtXrWeUtm9iZJu0spf03S\ne7X6dHoA48S8ADALZgWAWTArgIFyvYlRSrlf0p91fMlVWn1LoEopD0l6mZm9wjMTgJiYFwBmwawA\nMAtmBTBcrZ+Jca6kr6/7/KnpzwHAZswLALNgVgCYBbMCSCrNd4wwM54oCzgppQzqW+UyLwAfzAoA\ns2BWAJhF31nR+p0YT2n1yfSn/MD057ZUSknx8a53vWvpv6eUJ+uye+3bTbY10PcjkcHNC2ZFrF6Z\nFS/UTxqpZsXJkye15y17pBsl3STpRmnPW/bo5MmTTdfKPPthUVln7aZFzmXPi81ZvbuZb+2kwaxY\n2DVfzLr27maej8yvLaLOi3ks4yaGTT+2cq+kd0qSmV0i6dullD9aQiYAMTEvAMxiMLPirs/cpWM7\nj63925h07KXHdPdn726aKwK6qaObmTErRoBuug2xH9ebGGb2m5L+P0k/bGZ/aGbvNrP3mtl7JKmU\n8luSjpvZ1yT9uqT/2TPPsuzatat1hJmR1UemrFGMcV5kWidk9ZEpaxRDmxUPfPEBXfTdi7T3+F7p\nVmnv8b266ORFuv/h+zd8Xaa1sqiss3bTV+ZOvbsZAmZFfFlmhZS71yHOC5v3rRzLYmYlS9bJZKKV\nlZWl/p5mUp96WmTtq2/Wvt3MI1OvZqYywLOrGeYFs8IHs8IHs8LftOMtf23Za2We/eCRtaubvubJ\nuex50ZXVo5t5MCv8RZoVq3nivLbw2g9DeW0RaV7MMytaPxMDAAAAAABgJtzEAAAAAAAAKXCcZCBa\nvFUpC7rpxts+x4X9UEc33ZgV/mK9zTfWfojUjRSrn3jdMCu8xbvm7IeaSN1IsfrhOAkAAAAAABg8\nbmI4mEwmrSPMjKw+MmVFO5nWCVl9ZMqKtjKtlSxZs+SUcmVFW5nWCll9ZMraFzcxAAAAAABACjwT\nYyCinbeKhG66cXZ1XNgPdXTTjVnhL9ZZ5Vj7IVI3Uqx+4nXDrPAW75qzH2oidSPF6odnYgAAAAAA\ngMHjJoaDTOeQyOojU1a0k2mdkNVHpqxoK9NayZI1S04pV1a0lWmtkNVHpqx9cRMDAAAAAACkwDMx\nBiLaeatI6KYbZ1fHhf1QRzfdmBX+Yp1VjrUfInUjxeonXjfMCm/xrjn7oSZSN1KsfngmBgAAAAAA\nGDxuYjjIdA6JrD4yZUU7mdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEAAAAAAEAKPBNj\nIKKdt4qEbrpxdnVc2A91dNONWeEv1lnlWPshUjdSrH7idcOs8BbvmrMfaiJ1I8Xqh2diAAAAAACA\nweMmhoNM55DI6iNTVrSTaZ2Q1UemrGgr01rJkjVLTilXVrSVaa2Q1UemrH1xEwMAAAAAAKTAMzEG\nItp5q0jophtnV8eF/VBHN92YFf5inVWOtR8idSPF6ideN8wKb/GuOfuhJlI3Uqx+eCYGAAAAAAAY\nPG5iOMh0DomsPjJlRTuZ1glZfWTKirYyrZUsWbPklHJlRVuZ1gpZfWTK2hc3MQAAAAAAQAo8E2Mg\nop23ioRuunF2dVzYD3V0041Z4S/WWeVY+yFSN1KsfuJ1w6zwFu+asx9qInUjxeqHZ2IAAAAAAIDB\n4yaGg0znkMjqI1NWtJNpnZDVR6asaCvTWsmSNUtOKVdWtJVprZDVR6asfXETAwAAAAAApMAzMQYi\n2nmrSOimG2dXx4X9UEc33ZgV/mKdVY61HyJ1I8XqJ143zApv8a45+6EmUjdSrH54JgYAAAAAABg8\nbmI4yHQOiaw+MmVFO5nWCVl9ZMqKtjKtlSxZs+SUcmVFW5nWCll9ZMraFzcxAAAAAABACjwTYyCi\nnbeKhG66cXZ1XNgPdXTTjVnhL9ZZ5Vj7IVI3Uqx+4nXDrPAW75qzH2oidSPF6odnYgAAAAAAgMHj\nJoaDTOeQyOojU1a0k2mdkNVHpqxoK9NayZI1S04pV1a0lWmtkNVHpqx9cRMDAAAAAACkwDMxBiLa\neatI6KYbZ1fHhf1QRzfdmBX+Yp1VjrUfInUjxeonXjfMCm/xrjn7oSZSN1KsfngmBgAAAAAAGDxu\nYjjIdA6JrD4yZUU7mdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEAAAAAAEAKPBNjIKKd\nt4qEbrpxdnVc2A91dNONWeEv1lnlWPshUjdSrH7idcOs8BbvmrMfaiJ1I8Xqh2diLNGTx4/r0IED\nunHfPh06cEBPHj/eLMv7b3i/Dt97OMxCpJtukfqBv0jXm/1QRzfdIvYzRKeu+WUS+2GTKN1I8fqJ\n1A2WI9I1Zz/URetGitPPwroppaT4kFROnjxZWjrxxBPlA7t3l6dXb6iVp6Xygd27y4knntjwdUeO\nHFlKnsvfeXk58+CZZc/Ve4pUevWzqKyzdjOP7WRdRDfz2Jx1Gf30tToG2u/xRX5M/52aYVZ0894P\nzIrF2djPMGcFry3WLGo/ZHltsd2cLedF9Fmx3lBnRWuRZkUpsV5b8PeQjSLPi0W9rkj1ToxL33pp\n07tat11/vQ49/rjOmn5+lqRDjz+u266/vkkeM9Ozu57VQ699SFLbfuimW7R+4Cva9WY/1NFNt839\nDBHXfA37oVukfqJ1A3/Rrjn7oS5SN1Ksfhb2uqLv3Y9lf0gqulFFb1DRi1Wk5X9cNr17tfnjhn37\n5rwn1c/ed+0tukmrHypN+6GbXP2sJ/6LycLdsLKy5fW+rMHak1R0np63H3ZevrNce9O19EM32+hn\nmLOCPx/W8Gdnnn6idbPeUGdFa7H/fChN//ykm0T9LOh1Rap3Yux5ZI8OHzqsk8+ebDLArti/X89s\nyvSMpB3nnNOijjVl9X9a9kM3SfuBix3nnrvl9b5i//4m62/vyt7VEOv2w60fvFU333DzMms5LVI/\ndLO9foaIPx+2wJ+d3QL0E7YbuIn+50PLPz/pplukfhb2uqLFhe3zIXFudbPTZ4rewlm0zRbRzTwi\nn0XbTHPcBY36Mf13aoZZ0S3SOXdmRbeN/QxzVvDaYs2i9kOW1xa9n4nRYF5EnxXrDXVWtBZpVpQS\n67UFfw/ZKPK8WNTrilTvxDBr+92azjv/fF1z3326Zf9+/YSkW/bv1zX33afzzj+/SZ4Lz7tQt7/5\ndj346Qclte2HbrpF6we+ol1v9kMd3XTb3M8QRbrmN+7bx35YJ1I3Uqx+onUDf9GuOfuhLlI3Uqx+\nFvW6wkrJ8R5Rvj9zt0jfg5huusXrh+/n7ine9WY/1NBNN2bFuETbD9HQTx2zYnzYD3V0UzfPrEj1\nTgwAAAAAADBe3MRwMJlMWkeYGVl9ZMqKdjKtE7L6yJQVbWVaK1myZskp5cqKtjKtFbL6yJS1L25i\nAAAAAACAFHgmRk/xzirHOW9FN93i9cPZVU/xrjf7oYZuujErxiXafoiGfuqYFePDfqijmzqeiQEA\nAAAAAAaPmxgOMp1DIquPTFnRTqZ1QlYfmbKirUxrJUvWLDmlXFnRVqa1QlYfmbL2xU0MAAAAAACQ\nAs/E6CneWeU4563oplu8fji76ine9WY/1NBNN2bFuETbD9HQTx2zYnzYD3V0U8czMQAAAAAAwOBx\nE8NBpnNIZPWRKSvaybROyOojU1a0lWmtZMmaJaeUKyvayrRWyOojU9a+uIkBAAAAAABS4JkYPcU7\nqxznvBXddIvXD2dXPcW73uyHGrrpxqwYl2j7IRr6qWNWjA/7oY5u6ngmBgAAAAAAGDxuYjjIdA6J\nrD4yZUU7mdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEAAAAAAEAKPBOjp3hnleOct6Kb\nbvH64eyqp3jXm/1QQzfdmBXjEm0/REM/dcyK8WE/1NFNHc/EAAAAAAAAg8dNDAeZziGR1UemrGgn\n0zohq49MWdFWprWSJWuWnFKurGgr01ohq49MWftyv4lhZlea2VfN7FEzu3aLX/8+M7vXzI6a2TEz\nO+idCUA8zAoAs2BWAJgV8wIYJtdnYpjZDkmPSvqbkr4p6WFJbyulfHXd13xI0veVUj5kZi+X9B8k\nvaKU8peb/lmhzqLFO6sc57wV3XSL10/7s6uLnBXTrw0zL+Jdb/ZDDd10Y1aMS7T9EA391EWYFdMc\ng/17SDTshzq6qYv8TIyLJT1WSnmylPKcpDskXbXpa4qkndMf75T0p1u90AAwaMwKALNgVgCYFfMC\nGCjvmxjnSvr6us+/Mf259T4m6dVm9k1JX5b0i86Z3GU6h0RWH5myBsGsCI6sPjJlDWKUs0LKtVay\nZM2SU8qVNZBRzotMa4WsPjJl7euM1gEkvVHSl0opP2lmuyXdZ2avK6U8vfkLDx48qF27dkmSzj77\nbF1wwQVaWVmRtHaxlvX5qZ9r9fs/P89Ek8n2///r/11a5p/l86NHj4bK0/X50aNHQ+VZ//lkMtFt\nt90mSaf3UxIzzwop1rxg/y338777T4qRP8r6OPXjEydOKJnBzYpsn58SJc+iZwWfb/z81I8Tzgop\n6d9Don0uTTSZDHdWTEbyOijTrPB+JsYlkm4qpVw5/fw6SaWU8uF1X/NZSTeXUh6Yfv6vJV1bSvni\npn9WqLNo8c4qxzlvRTfd4vXT/uzqImfF9NfCzIt415v9UEM33ZgV4xJtP0RDP3URZsU0x2D/HhIN\n+6GObuoiPxPjYUmvNLPzzOxFkt4m6d5NX/OkpCskycxeIemHJT3hnAtALMwKALNgVgCYFfMCGCjX\nmxillO9Kep+kz0l6RNIdpZSvmNl7zew90y/7ZUk/bma/J+k+SX+/lPKfPHN52/wWqcjI6iNT1giY\nFfGR1UemrBGMdVZIudZKlqxZckq5skYx1nmRaa2Q1UemrH25PxOjlPI7kn5k08/9+roff0ur59EA\njBizAsAsmBUAZsW8AIbJ9ZkYixTtLFq8s8pxzlvRTbd4/cQ4u7pIkeZFvOvNfqihm27MinGJth+i\noZ86ZsX4sB/q6KYu8jMxAAAAAAAAFoKbGA4ynUMiq49MWdFOpnVCVh+ZsqKtTGslS9YsOaVcWdFW\nprVCVh+ZsvbFTQwAAAAAAJACz8ToKd5Z5TjnreimW7x+OLvqKd71Zj/U0E03ZsW4RNsP0dBPHbNi\nfNgPdXRTxzMxAAAAAADA4HETw0Gmc0hk9ZEpK9rJtE7I6iNTVrSVaa1kyZolp5QrK9rKtFbI6iNT\n1r64iQEAAAAAAFLgmRg9xTurHOe8Fd10i9cPZ1c9xbve7IcauunGrBiXaPshGvqpY1aMD/uhjm7q\neCYGAAAAAAAYPG5iOMh0DomsPjJlRTuZ1glZfWTKirYyrZUsWbPklHJlRVuZ1gpZfWTK2hc3MXoo\npUgvVqi3+UZBN93oZ1y43t3op45uxqeUousOXcc13wLd1NHN+HDN6+im26D6KaWk+FiNGsOd99xZ\n9AaVw/cebh3ltCj10E23mP2olAB7fJEfUeZFzOvdOsGaaP3QTTdmha8777mz7Lx8Z5hrHqiacN2U\nEqefmN0wKzzFvOatE6yim27R+plnVvBgz20qpejSn7tUD73mIe15ZI8e/PSDMmv/7KIID42hm25x\n++EBXB7iXm/2Qw3ddGNW+Il4zdkP3SL0E7cbZoWXuNec/VAToRspZj/zzIpUNzGkHFmliaSVxhlm\nNRFZPUyUJ+swX2zkmBcT5VknE5HVw0R5sjIr2pooz1qZKEfWiXLklHJlZVa0NVGetTIRWT1MlCNr\n/1lxxqKjeGp9v2X9HSyZpKIt72RNJtLKynKz9b3Lt6iss3Yzj75ZW9wB3Zx1Gf30FeAmtYuW84JZ\n0c17PzArfDArvH7/mPNinv2Q5bXFPDmXPS+YFe0xK+pav7bg7yEbZZkX8/zWPNhzG+76zF06tvPY\n6sWXJJOOvfSY7v7s3Ru+bmXZk2MOi8o6azfzyNzrMvpBHMyKbt77IXOvzIrxYV7UMSvWMCvArKjj\n7yEbjWFepHonRmsPfPEBXfTdi2TH124blVJ0/8P36+qfubphsvbophv9jAvXuxv91NHN+Ky/5p+f\nfF57V/Zyzafopo5uxodrXkc33YbYT6pnYmTJOplMln63rv/buJafta++Wdu8jStPrzyAqx1mhQ9m\nhQ9mhb9px1v+2rLXynzHSRaftaubvubJufzjJPWsHt3Mg1nhL9KsWM0T57WF134YymuLSPNinlnB\ncRIAAAAAAJAC78QYiCjfviciuunGfzEZF/ZDHd10Y1b4i/VfyGLth0jdSLH6idcNs8JbvGvOfqiJ\n1I0Uqx/eiQEAAAAAAAaPmxgOJpNJ6wgzI6uPTFnRTqZ1QlYfmbKirUxrJUvWLDmlXFnRVqa1QlYf\nmbL2xU0MAAAAAACQAs/EGIho560ioZtunF0dF/ZDHd10Y1b4i3VWOdZ+iNSNFKufeN0wK7zFu+bs\nh5pI3Uix+uGZGAAAAAAAYPC4ieEg0zkksvrIlBXtZFonZPWRKSvayrRWsmTNklPKlRVtZVorZPWR\nKWtf3MQAAAAAAAAp8EyMgYh23ioSuunG2dVxYT/U0U03ZoW/WGeVY+2HSN1IsfqJ1w2zwlu8a85+\nqInUjRSrH56JAQAAAAAABo+bGA4ynUMiq49MWdFOpnVCVh+ZsqKtTGslS9YsOaVcWdFWprVCVh+Z\nsvbFTQwAAAAAAJACz8QYiGjnrSKhm26cXR0X9kMd3XRjVviLdVY51n6I1I0Uq5943TArvMW75uyH\nmkjdSLH64ZkYAAAAAABg8LiJ4SDTOSSy+siUFe1kWidk9ZEpK9rKtFayZM2SU8qVFW1lWitk9ZEp\na1/cxAAAAAAAACnwTIyBiHbeKhK66cbZ1XFhP9TRTTdmhb9YZ5Vj7YdI3Uix+onXDbPCW7xrzn6o\nidSNFKsfnokBAAAAAAAGj5sYDjKdQyKrj0xZ0U6mdUJWH5myoq1MayVL1iw5pVxZ0VamtUJWH5my\n9sVNDAAAAAAAkALPxBiIaOetIqGbbpxdHRf2Qx3ddGNW+It1VjnWfojUjRSrn3jdMCu8xbvm7Iea\nSN1IsfrhmRgAAAAAAGDwuInhINM5JLL6yJQV7WRaJ2T1kSkr2sq0VrJkzZJTypUVbWVaK2T1kSlr\nX9zEAAAAAAAAKfBMjIGIdt4qErrpxtnVcWE/1NFNN2aFv1hnlWPth0jdSLH6idcNs8JbvGvOfqiJ\n1I0Uqx+eiQEAAAAAAAaPmxgOMp1DIquPTFnRTqZ1QlYfmbKirUxrJUvWLDmlXFnRVqa1QlYfmbL2\nxU0MAAAAAACQAs/EGIho560ioZtunF0dF/ZDHd10Y1b4i3VWOdZ+iNSNFKufeN0wK7zFu+bsh5pI\n3Uix+uGZGAAAAAAAYPC4ieEg0zkksvrIlBXtZFonZPWRKSvayrRWsmTNklPKlRVtZVorZPWRKWtf\n3MQAAAAAAAAp8EyMgYh23ioSuunG2dVxYT/U0U03ZoW/WGeVY+2HSN1IsfqJ1w2zwlu8a85+qInU\njRSrH56JAQAAAAAABo+bGA4ynUMiq49MWdFOpnVCVh+ZsqKtTGslS9YsOaVcWdFWprVCVh+ZsvbF\nTQwAAAAAAJACz8QYiGjnrSKhm26cXR0X9kMd3XRjVviLdVY51n6I1I0Uq5943TArvMW75uyHmkjd\nSLH64ZkYS/Tk8eM6dOCAbty3T4cOHNCTx483y/L+G96vw/ceDrMQ6aZbpH7gL9L1Zj/U0U23iP0M\n0alrfpnEftgkSjdSvH4idYPliHTN2Q910bqR4vSzsG5KKSk+JJWTJ0+Wlk488UT5wO7d5enVG2rl\naal8YPfucuKJJzZ83ZEjR5aS5/J3Xl7OPHhm2XP1niKVXv0sKuus3cxjO1kX0c08NmddRj99rY6B\n9nt8kR/Tf6dmmBXdvPcDs2JxNvYzzFnBa4s1i9oPWV5bbDdny3kRfVasN9RZ0VqkWVFKrNcW/D1k\no8jzYlGvK1K9E+PSt17a9K7Wbddfr0OPP66zpp+fJenQ44/rtuuvb5LHzPTsrmf10GsfktS2H7rp\nFq0f+Ip2vdkPdXTTbXM/Q8Q1X8N+6Bapn2jdwF+0a85+qIvUjRSrn4W9ruh792PZH5KKblTRG1T0\nYhVp+R+XTe9ebf64Yd++Oe9J9bP3XXuLbtLqh0rTfugmVz/rif9isnA3rKxseb0va7D2JBWdp+ft\nh52X7yzX3nQt/dDNNvoZ5qzgz4c1/NmZp59o3aw31FnRWuw/H0rTPz/pJlE/C3pdkeqdGHse2aPD\nhw7r5LMnmwywK/bv1zObMj0jacc557SoY01Z/Z+W/dBN0n7gYse55255va/Yv7/J+tu7snc1xLr9\ncOsHb9XNN9y8zFpOi9QP3WyvnyHiz4ct8GdntwD9hO0GbqL/+dDyz0+66Rapn4W9rmhxYft8SJxb\n3ez0maK3cBZts0V0M4/IZ9E20xx3QaN+TP+dmmFWdIt0zp1Z0W1jP8OcFby2WLOo/ZDltUXvZ2I0\nmBfRZ8V6Q50VrUWaFaXEem3B30M2ijwvFvW6ItU7Mczafrem884/X9fcd59u2b9fPyHplv37dc19\n9+m8889vkufC8y7U7W++XQ9++kFJbfuhm27R+oGvaNeb/VBHN9029zNEka75jfv2sR/WidSNFKuf\naN3AX7Rrzn6oi9SNFKufRb2usFJyvEeU78/cLdL3IKabbvH64fu5e4p3vdkPNXTTjVkxLtH2QzT0\nU8esGB/2Qx3d1M0zK9zfiWFmV5rZV83sUTO7tvI1K2b2JTP7fTM74p0JQDzMCgCzYFYAmBXzAhgm\n15sYZrZD0sckvVHSayS93cxetelrXibpVyX9dCnltZLe6plpGSaTSesIMyOrj0xZI2BWxEdWH5my\nRjDWWSHlWitZsmbJKeXKGsVY50WmtUJWH5my9uX9ToyLJT1WSnmylPKcpDskXbXpa35e0l2llKck\nqZTyJ86ZAMTDrAAwC2YFgFkxL4CBcn0mhpldLemNpZT3TD8/IOniUsrfXfc1/1TS92r1DulLJX20\nlPLJLf5Zoc6ixTurHOe8Fd10i9dP+7Ori5wV068NMy/iXW/2Qw3ddGNWjEu0/RAN/dRFmBXTHIP9\ne0g07Ic6uqmbZ1acsegwPZwh6UJJPynpLEkPmtmDpZSvbf7CgwcPateuXZKks88+WxdccIFWVlYk\nrb1tZlmfn/q5Vr//8/NMNJm0+/35PM/nk8lEt912mySd3k9JzDwrpFjzgs9zfC5Fy6Pmv/9kMtGJ\nEyeUDLOCz/l8iZ+f+nHCWSEl/XtItM+liSaTOHn4PObnp368kFnR93uzzvIh6RJJv7Pu8+skXbvp\na66VdOO6z/9PSVdv8c8qkXTlWdb3Z16vbz0eWb2uVd+sLZZOV9aga7n1919f2KwoweYFs6Kbx7Vi\nVvhgVrS17HkxTzUtZlsf8+Rc9tLJ0mkpMWZFWfC8YFZ0i/bawsNQXltEMs+s2DHvTZAX8LCkV5rZ\neWb2Iklvk3Tvpq+5R9JlZvY9ZnampD2SvuKcC0AszAoAs2BWAJgV8wIYKNdnYkir39pI0q9o9SGi\nnyil/GMze69W77x8fPo1H5T0bknflfQbpZR/vsU/p3hn3Y54Z5XjnLeim27x+glzdnUhs2L6dWHm\nRbzrzX7gV66OAAAgAElEQVSooZtuzIpxibYfoqGfuiizQhru30OiYT/U0U3dPLPC/SbGokQbHvFe\nXMbZIHTTLV4/cV5sLEqkeRHverMfauimG7NiXKLth2jop45ZMT7shzq6qZtnVngfJxml9Q8viY6s\nPjJlRTuZ1glZfWTKirYyrZUsWbPklHJlRVuZ1gpZfWTK2hc3MQAAAAAAQAocJ+kp3tt847xViW66\nxeuHt316ine92Q81dNONWTEu0fZDNPRTx6wYH/ZDHd3UcZwEAAAAAAAMHjcxHGQ6h0RWH5myop1M\n64SsPjJlRVuZ1kqWrFlySrmyoq1Ma4WsPjJl7YubGAAAAAAAIAWeidFTvLPKcc5b0U23eP1wdtVT\nvOvNfqihm27MinGJth+ioZ86ZsX4sB/q6KaOZ2IAAAAAAIDB4yaGg0znkMjqI1NWtJNpnZDVR6as\naCvTWsmSNUtOKVdWtJVprZDVR6asfXETAwAAAAAApMAzMXqKd1Y5znkruukWrx/OrnqKd73ZDzV0\n041ZMS7R9kM09FPHrBgf9kMd3dTxTAwAAAAAADB43MRwkOkcEll9ZMqKdjKtE7L6yJQVbWVaK1my\nZskp5cqKtjKtFbL6yJS1L25iAAAAAACAFHgmRk/xzirHOW9FN93i9cPZVU/xrjf7oYZuujErxiXa\nfoiGfuqYFePDfqijmzqeiQEAAAAAAAaPmxgOMp1DIquPTFnRTqZ1QlYfmbKirUxrJUvWLDmlXFnR\nVqa1QlYfmbL2xU0MAAAAAACQAs/E6CneWeU4563oplu8fji76ine9WY/1NBNN2bFuETbD9HQTx2z\nYnzYD3V0U8czMQAAAAAAwOBxE8NBpnNIZPWRKSvaybROyOojU1a0lWmtZMmaJaeUKyvayrRWyOoj\nU9a+uIkBAAAAAABS4JkYPcU7qxznvBXddIvXD2dXPcW73uyHGrrpxqwYl2j7IRr6qWNWjA/7oY5u\n6ngmBgAAAAAAGDxuYjjIdA6JrD4yZUU7mdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEA\nAAAAAEAKPBOjp3hnleOct6KbbvH64eyqp3jXm/1QQzfdmBXjEm0/REM/dcyK8WE/1NFNHc/EAAAA\nAAAAg8dNDAeZziGR1UemrGgn0zohq49MWdFWprWSJWuWnFKurGgr01ohq49MWfvqvIlhZt9nZru3\n+PnX+UWKr5QivVih3uYbBd10G3I/zIvnG/L1XgT6qRtyN8yKrZVSdN2h6wZ5zedFN3VD7oZZsbUh\nX/N50U23QfVTStnyQ9LPSfqmpKOSHpH0Y+t+7d/X/n9eH6tRY7jznjuL3qBy+N7DraOcFqUeuukW\nsx+VMv/+ZF5sIeb1bp1gTbR+6KYbs8LXnffcWXZevjPMNQ9UTbhuSonTT8xumBWeYl7z1glW0U23\naP3MMyu63onxv0r6G6WUCyS9W9InzexvT39tUA/r2Y5Sim755C3S35I+cvtHTg02iG5eyMD7YV5s\nMvDrPTf6qRt4N8yKLZy65t/Z950hXvO50E3dwLthVmxh4Nd8LnTTbWj9VL87iZkdK6X86LrP/6qk\nz0r6vyQdLKVcuJyIp3//ImUpeyJppXGGWU1EVg8T5ck6/1PEmRd9TZRnnUxEVg8T5cnKrGhrojxr\nZaIcWSfKkVPKlZVZ0dZEedbKRGT1MFGOrP1nxRkdv/YdM9tdSnlckkop3zKzFUn/StJr+vxm82p9\nw6iUokt/7lI99JqHVu8BF2nPI3v04KcflNla/5OJtLKy3Gx9v33PorLO2s08+mZt8a2NNmddRj99\nLei3Z15s+L2ZFV289wOzwgezwuv3jzkv5tkPWV5bzJNz2fOCWcGsiDorpPavLfh7yEZZ5sU8v3XX\ncZL/SdIOM3v1qZ8opXxH0pWS/k7/3zKvuz5zl47tPLb2JjaTjr30mO7+7N0bvm5l2ZNjDovKOms3\n88jc6zL6aYx5sQ6zopv3fsjcK7NifJgXdcyKNcwKZgWzoo6/h2w0hnlRfSdGKeXLkmRmv29mn5T0\nTyS9ePq/F0n65FISBvLAFx/QRd+9SHZ87bZRKUX3P3y/rv6Zqxsma49uug29H+bFRkO/3vOin7qh\nd8OseL711/zzk89r78reQV3zedBN3dC7YVY839Cv+TzoptsQ+6k+E+P0F5idJenDkv6GpJ2SPiXp\nw6WUk/7xNuQoWR5AMplMln63rv/buJafta++Wdu8jStPr2bzn11d989iXmwDs8IHs8IHs8LftOMt\nf23Za2W+4ySLz9rVTV/z5Fz+cZJ6Vo9u5sGs8BdpVqzmifPawms/DOW1RaR5Mc+s6DpOcspzkv5C\n0ku0egf0+LIHB4A0mBcAZsGsADALZgWA55nlnRhflnSPpH8o6eWSfk3Sfy2lvNU/3oYcoe6ARtPi\nLl8WdNNtwf/FhHkRHPuhjm66MSv8xfovZLH2Q6RupFj9xOuGWeEt3jVnP9RE6kaK1c88s2KWmxgX\nlVK+uOnn3lFKWepZtGjDI5poGyQSuum24BcbzIvg2A91dNONWeEv1ovLWPshUjdSrH7idcOs8Bbv\nmrMfaiJ1I8Xqx/U4yebBMf250T1MZzsmk0nrCDMjq49MWReJebE9mdYJWX1kyrpIzIrty7RWsmTN\nklPKlXWRmBXbl2mtkNVHpqx9zfJMDAAAAAAAgOZe8DhJFNHexhVNtLcqRUI33Rb5ts8omBd17Ic6\nuunGrPAX622+sfZDpG6kWP3E64ZZ4S3eNWc/1ETqRorVj/d3JwEAAAAAAGiOmxgOMp1DIquPTFnR\nTqZ1QlYfmbKirUxrJUvWLDmlXFnRVqa1QlYfmbL2xU0MAAAAAACQAs/EGIho560ioZtunF0dF/ZD\nHd10Y1b4i3VWOdZ+iNSNFKufeN0wK7zFu+bsh5pI3Uix+uGZGAAAAAAAYPC4ieEg0zkksvrIlBXt\nZFonZPWRKSvayrRWsmTNklPKlRVtZVorZPWRKWtf3MQAAAAAAAAp8EyMgYh23ioSuunG2dVxYT/U\n0U03ZoW/WGeVY+2HSN1IsfqJ1w2zwlu8a85+qInUjRSrH56JAQAAAAAABo+bGA4ynUMiq49MWdFO\npnVCVh+ZsqKtTGslS9YsOaVcWdFWprVCVh+ZsvbFTQwAAAAAAJACz8QYiGjnrSKhm26cXR0X9kMd\n3XRjVviLdVY51n6I1I0Uq5943TArvMW75uyHmkjdSLH64ZkYAAAAAABg8LiJ4SDTOSSy+siUFe1k\nWidk9ZEpK9rKtFayZM2SU8qVFW1lWitk9ZEpa1/cxAAAAAAAACnwTIyBiHbeKhK66cbZ1XFhP9TR\nTTdmhb9YZ5Vj7YdI3Uix+onXDbPCW7xrzn6oidSNFKsfnokBAAAAAAAGj5sYDjKdQyKrj0xZ0U6m\ndUJWH5myoq1MayVL1iw5pVxZ0VamtUJWH5my9sVNDAAAAAAAkALPxBiIaOetIqGbbpxdHRf2Qx3d\ndGNW+It1VjnWfojUjRSrn3jdMCu8xbvm7IeaSN1IsfrhmRgAAAAAAGDwuInhINM5JLL6yJQV7WRa\nJ2T1kSkr2sq0VrJkzZJTypUVbWVaK2T1kSlrX+43MczsSjP7qpk9ambXdnzdj5nZc2b2Zu9MAOJh\nVgCYBbMCwKyYF8AwuT4Tw8x2SHpU0t+U9E1JD0t6Wynlq1t83X2S/kLSvyyl3L3FPyvUWbRoop23\nioRuukU4u7rIWTH9OuZFBfuhjm66MSv8xTqrHGs/ROpGitVPvG7az4ppjsH+PSTeNWc/1ETqRorV\nT+RnYlws6bFSypOllOck3SHpqi2+7hpJhyX9sXMeADExKwDMglkBYFbMC2CgvG9inCvp6+s+/8b0\n504zs3Mk/Wwp5V9Ian7XdhEynUMiq49MWYNgVgRHVh+ZsgYxylkh5VorWbJmySnlyhrIKOdFprVC\nVh+ZsvZ1RusAkv6ZpPVn1KoD5ODBg9q1a5ck6eyzz9YFF1yglZUVSWsXa6yfSxNNJtv//5/SOv8s\nnx89ejRUnq7Pjx49GirP+s8nk4luu+02STq9n5KYeVZIOebFKVHyjH3/STHyR/n81I9PnDihZNLO\nilM/1/raj2W2Rf6zejufn8Ks6CXl30NO/VzrtbeWZ6LJJN/a3c7nY3gdlGlWeD8T4xJJN5VSrpx+\nfp2kUkr58LqveeLUDyW9XNIzkt5TSrl30z8r1Fm0aKKdt4qEbrpFOLu6yFkx/VrmRQX7oY5uujEr\n/MU6qxxrP0TqRorVT7xu2s+KaY7B/j0k3jVnP9RE6kaK1c88s8L7Jsb3SPoPWn2gzrck/TtJby+l\nfKXy9bdK+kyGB+pEE22DREI33SK82FjkrJj+OvOigv1QRzfdmBX+Yr24jLUfInUjxeonXjftZ8U0\nx2D/HhLvmrMfaiJ1I8XqJ+yDPUsp35X0Pkmfk/SIpDtKKV8xs/ea2Xu2+r945lmWzW+RioysPjJl\njYBZER9ZfWTKGsFYZ4WUa61kyZolp5QraxRjnReZ1gpZfWTK2pf7MzFKKb8j6Uc2/dyvV772F7zz\nAIiJWQFgFswKALNiXgDD5HqcZJGivY0rmmhvVYqEbrpFedvnIjEv6tgPdXTTjVnhL9bbfGPth0jd\nSLH6idcNs8JbvGvOfqiJ1I0Uq5+wx0mG6Mnjx3XowAHduG+fDh04oCePH2+W5f03vF+H7z0cZiHS\nTbdI/cBfpOvNfqijm24R+xmiU9f8Mon9sEmUbqR4/UTqBssR6ZqzH+qidSPF6Wdh3ZRSUnxIKidP\nniwtnXjiifKB3bvL06s31MrTUvnA7t3lxBNPbPi6I0eOLCXP5e+8vJx58Myy5+o9RSq9+llU1lm7\nmcd2si6im3lszrqMfvpaHQPt9/giP6b/Ts0wK7p57wdmxeJs7GeYs4LXFmsWtR+yvLbYbs6W8yL6\nrFhvqLOitUizopRYry34e8hGkefFol5XpHonxqVvvbTpXa3brr9ehx5/XGdNPz9L0qHHH9dt11/f\nJI+Z6dldz+qh1z4kqW0/dNMtWj/wFe16sx/q6Kbb5n6GiGu+hv3QLVI/0bqBv2jXnP1QF6kbKVY/\nC3td0ffux7I/JBXdqKI3qOjFKtLyPy6b3r3a/HHDvn1z3pPqZ++79hbdpNUPlab90E2uftYT/8Vk\n4W5YWdnyel/WYO1JKjpPz9sPOy/fWa696Vr6oZtt9DPMWcGfD2v4szNPP9G6WW+os6K12H8+lKZ/\nftJNon4W9Loi1Tsx9jyyR4cPHdbJZ082GWBX7N+vZzZlekbSjnPOaVHHmrL6Py37oZuk/cDFjnPP\n3fJ6X7F/f5P1t3dl72qIdfvh1g/eqptvuHmZtZwWqR+62V4/Q8SfD1vgz85uAfoJ2w3cRP/zoeWf\nn3TTLVI/C3td0eLC9vmQOLe62ekzRW/hLNpmi+hmHpHPom2mOe6CRv2Y/js1w6zoFumcO7Oi28Z+\nhjkreG2xZlH7Ictri97PxGgwL6LPivWGOitaizQrSon12oK/h2wUeV4s6nVFqndimLX9bk3nnX++\nrrnvPt2yf79+QtIt+/frmvvu03nnn98kz4XnXajb33y7Hvz0g5La9kM33aL1A1/Rrjf7oY5uum3u\nZ4giXfMb9+1jP6wTqRspVj/RuoG/aNec/VAXqRspVj+Lel1hpeR4jyjfn7lbpO9BTDfd4vXD93P3\nFO96sx9q6KYbs2Jcou2HaOinjlkxPuyHOrqpm2dWpHonBgAAAAAAGC9uYjiYTCatI8yMrD4yZUU7\nmdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEAAAAAAEAKPBOjp3hnleOct6KbbvH64eyq\np3jXm/1QQzfdmBXjEm0/REM/dcyK8WE/1NFNHc/EAAAAAAAAg8dNDAeZziGR1UemrGgn0zohq49M\nWdFWprWSJWuWnFKurGgr01ohq49MWfviJgYAAAAAAEiBZ2L0FO+scpzzVnTTLV4/nF31FO96sx9q\n6KYbs2Jcou2HaOinjlkxPuyHOrqp45kYAAAAAABg8LiJ4SDTOSSy+siUFe1kWidk9ZEpK9rKtFay\nZM2SU8qVFW1lWitk9ZEpa1/cxAAAAAAAACnwTIye4p1VjnPeim66xeuHs6ue4l1v9kMN3XRjVoxL\ntP0QDf3UMSvGh/1QRzd1PBMDAAAAAAAMHjcxHGQ6h0RWH5myop1M64SsPjJlRVuZ1kqWrFlySrmy\noq1Ma4WsPjJl7YubGAAAAAAAIAWeidFTvLPKcc5b0U23eP1wdtVTvOvNfqihm27MinGJth+ioZ86\nZsX4sB/q6KaOZ2IAAAAAAIDB4yaGg0znkMjqI1NWtJNpnZDVR6asaCvTWsmSNUtOKVdWtJVprZDV\nR6asfXETAwAAAAAApMAzMXqKd1Y5znkruukWrx/OrnqKd73ZDzV0041ZMS7R9kM09FPHrBgf9kMd\n3dTxTAwAAAAAADB43MRwkOkcEll9ZMqKdjKtE7L6yJQVbWVaK1myZskp5cqKtjKtFbL6yJS1L25i\nAAAAAACAFHgmRk/xzirHOW9FN93i9cPZVU/xrjf7oYZuujErxiXafoiGfuqYFePDfqijmzqeiQEA\nAAAAAAaPmxgOMp1DIquPTFnRTqZ1QlYfmbKirUxrJUvWLDmlXFnRVqa1QlYfmbL2xU0MAAAAAACQ\nAs/E6CneWeU4563oplu8fji76ine9WY/1NBNN2bFuETbD9HQTx2zYnzYD3V0U8czMQAAAAAAwOBx\nE8NBpnNIZPWRKSvaybROyOojU1a0lWmtZMmaJaeUKyvayrRWyOojU9a+uIkBAAAAAABS4JkYPcU7\nqxznvBXddIvXD2dXPcW73uyHGrrpxqwYl2j7IRr6qWNWjA/7oY5u6ngmBgAAAAAAGDxuYjjIdA6J\nrD4yZUU7mdYJWX1kyoq2Mq2VLFmz5JRyZUVbmdYKWX1kytoXNzEAAAAAAEAKPBOjp3hnleOct6Kb\nbvH64eyqp3jXm/1QQzfdmBXjEm0/REM/dcyK8WE/1NFNHc/EAAAAAAAAg8dNDAeZziGR1UemrGgn\n0zohq49MWdFWprWSJWuWnFKurGgr01ohq49MWfviJkYPpRTpxQr1Nt8o6KYb/YwL17sb/dTRzfiU\nUnTdoeu45lugmzq6GR+ueR3ddBtUP6WUFB+rUWO48547i96gcvjew62jnBalHrrpFrMflRJgjy/y\nI8q8iHm9WydYE60fuunGrPB15z13lp2X7wxzzQNVE66bUuL0E7MbZoWnmNe8dYJVdNMtWj/zzAoe\n7LlNpRRd+nOX6qHXPKQ9j+zRg59+UGbtn10U4aExdNMtbj88gMtD3OvNfqihm27MCj8Rrzn7oVuE\nfuJ2w6zwEveasx9qInQjxexnnlmR6iaGlCOrNJG00jjDrCYiq4eJ8mQd5ouNHPNiojzrZCKyepgo\nT1ZmRVsT5VkrE+XIOlGOnFKurMyKtibKs1YmIquHiXJk7T8rzlh0FE+t77esv4Mlk1S05Z2syURa\nWVlutr53+RaVddZu5tE3a4s7oJuzLqOfvgLcpHbRcl4wK7p57wdmhQ9mhdfvH3NezLMfsry2mCfn\nsucFs6I9ZkVd69cW/D1koyzzYp7fmgd7bsNdn7lLx3YeW734kmTSsZce092fvXvD160se3LMYVFZ\nZ+1mHpl7XUY/iINZ0c17P2TulVkxPsyLOmbFGmYFmBV1/D1kozHMi1TvxGjtgS8+oIu+e5Hs+Npt\no1KK7n/4fl39M1c3TNYe3XSjn3Hhenejnzq6GZ/11/zzk89r78pervkU3dTRzfhwzevoptsQ+0n1\nTIwsWSeTydLv1vV/G9fys/bVN2ubt3Hl6ZUHcLXDrPDBrPDBrPA37XjLX1v2WpnvOMnis3Z109c8\nOZd/nKSe1aObeTAr/EWaFat54ry28NoPQ3ltEWlezDMrOE4CAAAAAABS4J0YAxHl2/dERDfd+C8m\n48J+qKObbswKf7H+C1ms/RCpGylWP/G6YVZ4i3fN2Q81kbqRYvXDOzEAAAAAAMDgcRPDwWQyaR1h\nZmT1kSkr2sm0TsjqI1NWtJVprWTJmiWnlCsr2sq0VsjqI1PWvriJAQAAAAAAUuCZGAMR7bxVJHTT\njbOr48J+qKObbswKf7HOKsfaD5G6kWL1E68bZoW3eNec/VATqRspVj88EwMAAAAAAAweNzEcZDqH\nRFYfmbKinUzrhKw+MmVFW5nWSpasWXJKubKirUxrhaw+MmXti5sYAAAAAAAgBZ6JMRDRzltFQjfd\nOLs6LuyHOrrpxqzwF+uscqz9EKkbKVY/8bphVniLd83ZDzWRupFi9RP6mRhmdqWZfdXMHjWza7f4\n9Z83sy9PP+43sx/1zgQgHmYFgFkwKwDMinkBDJPrTQwz2yHpY5LeKOk1kt5uZq/a9GVPSLq8lPJ6\nSb8s6Tc8My1DpnNIZPWRKWsEzIr4yOojU9YIxjorpFxrJUvWLDmlXFmjGOu8yLRWyOojU9a+vN+J\ncbGkx0opT5ZSnpN0h6Sr1n9BKeULpZQ/n376BUnnOmcCEA+zAsAsmBUAZsW8AAbK9ZkYZna1pDeW\nUt4z/fyApItLKX+38vUflPTDp75+06+FOosWTbTzVpHQTbcIZ1cXOSumv868qGA/1NFNN2aFv1hn\nlWPth0jdSLH6iddN+1kxzTHYv4fEu+bsh5pI3Uix+plnVpyx6DB9mdk+Se+WdFntaw4ePKhdu3ZJ\nks4++2xdcMEFWllZkbT2tpmxfi5NNJnEycPncT+fTCa67bbbJOn0fspkllkhMS/4fPufS7HytP78\n1I9PnDihjDLOilM/1/ra83muz09hVvSX7e8hp36u9dpbyzPRZNJ+L/B57M9P/Xghs6KU4vYh6RJJ\nv7Pu8+skXbvF171O0mOSdnf8s0oWR44cWfrv2beeFln76pu1xdLJ1Ot0b7nOghf6WOSsKInmBbPC\nB7PCB7PCX1eeZa+VearxyOpxrebJueyl05U16DpuOivKgudF0I63NPbXFl7XaiivLSKt5XlmxY75\nb4N0eljSK83sPDN7kaS3Sbp3/ReY2Q9JukvSO0opjzvnARATswLALJgVAGbFvAAGyvWZGNLqtzaS\n9CtafYjoJ0op/9jM3qvVOy8fN7PfkPRmSU9KMknPlVIu3uKfU7yzZhbtvFUkdNMt0NnVhcyK6T+L\neVHBfqijm27MCn+xzirH2g+RupFi9ROvmxizQhru30PiXXP2Q02kbqRY/cwzK9xvYixKtOERTbQN\nEgnddIv0YmNRmBd17Ic6uunGrPAX68VlrP0QqRspVj/xumFWeIt3zdkPNZG6kWL1M8+s8D5OMkqb\nH7QUGVl9ZMqKdjKtE7L6yJQVbWVaK1myZskp5cqKtjKtFbL6yJS1L25iAAAAAACAFDhOMhDR3qoU\nCd10422f48J+qKObbswKf7He5htrP0TqRorVT7xumBXe4l1z9kNNpG6kWP1wnAQAAAAAAAweNzEc\nZDqHRFYfmbKinUzrhKw+MmVFW5nWSpasWXJKubKirUxrhaw+MmXti5sYAAAAAAAgBZ6JMRDRzltF\nQjfdOLs6LuyHOrrpxqzwF+uscqz9EKkbKVY/8bphVniLd83ZDzWRupFi9cMzMQAAAAAAwOBxE8NB\npnNIZPWRKSvaybROyOojU1a0lWmtZMmaJaeUKyvayrRWyOojU9a+uIkBAAAAAABS4JkYAxHtvFUk\ndNONs6vjwn6oo5tuzAp/sc4qx9oPkbqRYvUTrxtmhbd415z9UBOpGylWPzwTAwAAAAAADB43MRxk\nOodEVh+ZsqKdTOuErD4yZUVbmdZKlqxZckq5sqKtTGuFrD4yZe2LmxgAAAAAACAFnokxENHOW0VC\nN904uzou7Ic6uunGrPAX66xyrP0QqRspVj/xumFWeIt3zdkPNZG6kWL1wzMxAAAAAADA4HETw0Gm\nc0hk9ZEpK9rJtE7I6iNTVrSVaa1kyZolp5QrK9rKtFbI6iNT1r64iQEAAAAAAFLgmRgDEe28VSR0\n042zq+PCfqijm27MCn+xzirH2g+RupFi9ROvG2aFt3jXnP1QE6kbKVY/PBMDAAAAAAAMHjcxtunJ\n48d16MAB3bhvnw4dOKAnjx9/3tcs6xzS+294vw7fe3iuu2mLzDpLN/PYTtZFdDOPrbJ694NYmBXd\nPPcDs2JxWvczFqeu+WVS03mxqOvt8dqiq5u+tpuz5X7omhUe3SCmKLNCivfawns/DOW1Ret5sbBu\nSikpPiSVkydPlpZOPPFE+cDu3eXp1XcFlael8oHdu8uJJ57Y8HVHjhxZSp7L33l5OfPgmWXP1XuK\nVHr1s6iss3Yzj+1kXUQ389icdRn99LU6Btrv8UV+TP+dmmFWdPPeD8yKxdnYzzBnBa8t1ixqP2R5\nbbHdnC3nRfRZsd5QZ0VrkWZFKbFeW/D3kI0iz4tFva5oPhRmDiqVPVfvKXfec2ezFxw37d9/+uKX\ndYvgpv37m+TZ+669RTep6EYVqTTth266RetnPV5sLF60681+qKObbhv7Geas4JqvYT90i9RPtG7W\nG+qsaC3aNWc/1EXqppRY/SzqdUWq4yQPvfYhvfWmt2rHmTtkZkv/+N1PfUpnbcp0lqST3/xmizrW\nTB+H0rIfuknaD1ycfOqpLa/3737qU03W3+cnn18NsW4//ML//gv60C99aJm1nBapH7rZXj9DxJ8P\nW+DPzm4B+gnbDdxE//Oh5Z+fdNMtUj+Lel2R6ibGnkf26PChwzr57Mkmd2Gv2L9fz2zK9IykHeec\ns+Hnlv69eadHivr0c+TIkaV2M49evc7RzTwfm3tdRj+IY8e55255va/Yv99l/73Qx96Vvash1u2H\nWz94q26+4eaZ/50WOddm7WdR+8+7m3ls7tW7m3nXzhDx2mILc/7ZmeW1Re9OG7y24HUFeG1Rt4w/\nO9wWro0AAA1zSURBVHlt4bNuelt28L4f4tzq85w+U/QWzqJttohu5hH5LNpm4m2fC8es6BbpnDuz\notvGfoY5K3htsWZR+yHLa4vez8RoMC+iz4r1hjorWos0K0qJ9dqCv4dsFHleLOp1RfOhMHPQAMOj\nlNVFcNP+/eWy6Tmiln9Y/L3r/145fO/hcvLkyRKhHrrpFqmf9Xix4SPS9WY/1NFNt439MCu8nLrm\nN+zbx37YJEo3pcTrJ1I36zEr/ES65uyHumjdlBKnn0W9rrBScrxH1MxKpKxmplh5Vm+tRUA33eL1\nYyqlDOrEe6R5Ee96sx9q6KYbs2Jcou2HaOinjlkxPuyHOrqpm2dWpHomRhZLfybGHMjqI1NWtJNp\nnZDVR6asaCvTWsmSNUtOKVdWtJVprZDVR6asfXETAwAAAAAApMBxkp7ivc03zluV6KZbvH5426en\neNeb/VBDN92YFeMSbT9EQz91zIrxYT/U0U0dx0kAAAAAAMDgcRPDQaZzSGT1kSkr2sm0TsjqI1NW\ntJVprWTJmiWnlCsr2sq0VsjqI1PWvriJAQAAAAAAUuCZGD3FO6sc57wV3XSL1w9nVz3Fu97shxq6\n6casGJdo+yEa+qljVowP+6GObup4JgYAAAAAABg8bmI4yHQOiaw+MmVFO5nWCVl9ZMqKtjKtlSxZ\ns+SUcmVFW5nWCll9ZMraFzcxAAAAAABACjwTo6d4Z5XjnLeim27x+uHsqqd415v9UEM33ZgV4xJt\nP0RDP3XMivFhP9TRTR3PxAAAAAAAAIPHTQwHmc4hkdVHpqxoJ9M6IauPTFnRVqa1kiVrlpxSrqxo\nK9NaIauPTFn74iYGAAAAAABIgWdi9BTvrHKc81Z00y1eP5xd9RTverMfauimG7NiXKLth2jop45Z\nMT7shzq6qeOZGAAAAAAAYPC4ieEg0zkksvrIlBXtZFonZPWRKSvayrRWsmTNklPKlRVtZVorZPWR\nKWtf3MQAAAAAAAAp8EyMnuKdVY5z3opuusXrh7OrnuJdb/ZDDd10Y1aMS7T9EA391DErxof9UEc3\ndTwTAwAAAAAADB43MRxkOodEVh+ZsqKdTOuErD4yZUVbmdZKlqxZckq5sqKtTGuFrD4yZe2LmxgA\nAAAAACAFnonRU7yzynHOW9FNt3j9cHbVU7zrzX6ooZtuzIpxibYfoqGfOmbF+LAf6uimjmdiAAAA\nAACAweMmhoNM55DI6iNTVrSTaZ2Q1UemrGgr01rJkjVLTilXVrSVaa2Q1UemrH1xEwMAAAAAAKTA\nMzF6indWOc55K7rpFq8fzq56ine92Q81dNONWTEu0fZDNPRTx6wYH/ZDHd3U8UwMAAAAAAAweNzE\ncJDpHBJZfWTKinYyrROy+siUFW1lWitZsmbJKeXKirYyrRWy+siUtS9uYgAAAAAAgBR4JkZP8c4q\nxzlvRTfd4vXD2VVP8a43+6GGbroxK8Yl2n6Ihn7qmBXjw36oo5s6nokBAAAAAAAGj5sYDjKdQyKr\nj0xZ0U6mdUJWH5myoq1MayVL1iw5pVxZ0VamtUJWH5my9sVNDAAAAAAAkALPxOgp3lnlOOet6KZb\nvH44u+op3vVmP9TQTTdmxbhE2w/R0E8ds2J82A91dFPHMzEAAAAAAMDgcRPDQaZzSGT1kSkr2sm0\nTsjqI1NWtJVprWTJmiWnlCsr2sq0VsjqI1PWvtxvYpjZlWb2VTN71MyurXzNR83sMTM7amYXeGfy\ndvTo0dYRZkZWH5myRsGsiI2sPjJljWKMs0LKtVayZM2SU8qVNZIxzotMa4WsPjJl7cv1JoaZ7ZD0\nMUlvlPQaSW83s1dt+po3SdpdSvlrkt4r6dc8My1CKUV6sapnlb/97W8vOVF/i876Qt3MYwi9evaT\nGbMiPo+sXvthCL0yK7Y21FkxiyGs62iy5JRyZY1irPMi01ohq49MWfvyfifGxZIeK6U8WUp5TtId\nkq7a9DVXSbpdkkopD0l6mZm9wjnXXO76zF3SX5fu/uzdraOEQzfd6KeKWTFC9FNHN1WDnBUAXDAv\ngIHyvolxrqSvr/v8G9Of6/qap7b4mjBKKbrlk7dIf0v6yO0f2fK/kp04cWL5wXpaZNZZuplH9l69\n+0mOWRHcorN67ofsvTIrOg1uVswq+7qOKEtOKVfWQEY5LzKtFbL6yJS1t1KK24ekqyV9fN3nByR9\ndNPXfEbSj6/7/HclXbjFP6vwwQcfPh+ec2DZs4J5wQcffh/MCj744GOWj9azYtHzonWffPAx1I++\n+/sM+XpK0g+t+/wHpj+3+Wt+8AW+RmVg328awAYLmxUS8wIYMGYFgFnx9xBgoLyPkzws6ZVmdp6Z\nvUjS2yTdu+lr7pX0Tkkys0skfbuU8kfOuQDEwqwAMAtmBYBZMS+AgXJ9J0Yp5btm9j5Jn9PqDZNP\nlFK+YmbvXf3l8vFSym+Z2U+Z2dckPSPp3Z6ZAMTDrAAwC2YFgFkxL4Dhsuk5LwAAAAAAgNC8j5Ns\nm5ldaWZfNbNHzezaytd81MweM7OjZnbBsjOuy9GZ1cx+3sy+PP2438x+tEXOaZYX7HX6dT9mZs+Z\n2ZuXmW/d7z/L9V8xsy+Z2e+b2ZFlZ1yX44Wu//eZ2b3TdXrMzA42iHkqyyfM7I/M7Pc6vibEvpoV\ns8JHllkxzcC8WDBmBbNiVswKH8yKtpgXy8+57uuYFduQZVZMsyx+XrR+cvCmJ//ukPQ1SedJ+l5J\nRyW9atPXvEnS/z398R5JXwic9RJJL5v++MrIWdd93b+W9FlJb46YU9LLJD0i6dzp5y+P2qmkD0m6\n+VROSX8q6YxGeS+TdIGk36v8eoh9teD+Q/w7MSua9sq82H5WZkXg/cescOuVWbH9rIOaFdvoP8S/\nV5Z5waxomjXErJj+/gufF9HeiXGxpMdKKU+WUp6TdIekqzZ9zVWSbpekUspDkl5mZq9YbkxJM2Qt\npXyhlPLn00+/oHbfd3qWXiXpGkmHJf3xMsOtM0vOn5d0VynlKUkqpfzJkjOeMkvWImnn9Mc7Jf1p\nKeUvl5hxLUgp90v6s44vibKvZsWs8JFlVkjMCxfMCmbFjJgVPpgVbTEvFo9Z4SPNrJB85kW0mxjn\nSvr6us+/oedvuM1f89QWX7MMs2Rd7+9I+m3XRHUvmNXMzpH0s6WUfyGp1beRmqXTH5b0/WZ2xMwe\nNrN3LC3dRrNk/ZikV5vZNyV9WdIvLilbH1H21ayYFT6yzAqJedFKlH01K2aFD2aFD2ZFW8yLxWNW\n+BjSrJB67CvX706CVWa2T6tPO76sdZYO/0zS+vNUUb8f9hmSLpT0k5LOkvSgmT1YSvla21hbeqOk\nL5VSftLMdku6z8xeV0p5unUwxMSsWDjmBQaJWbFwzAoMVoJ5wazwMehZEe0mxlOSfmjd5z8w/bnN\nX/ODL/A1yzBLVpnZ6yR9XNKVpZSut9F4miXrRZLuMDPT6rmpN5nZc6WUzd9P29MsOb8h6U9KKf9F\n0n8xs38r6fVaPRe2TLNkfbekmyWplPK4mR2X9CpJX1xKwu2Jsq9mxazwkWVWSMyLVqLsq1kxK3ww\nK3wwK9piXiwes8LHkGaF1GdfzfpAjmV8SPoerT2k5EVafUjJX9/0NT+ltQd/XKJ2D9SZJesPSXpM\n0iXRe9309beqzYM9Z+n0VZLum37tmZKOSXp10Ky/KunG6Y9fodW3SX1/w3WwS9Kxyq+F2FcL7j/E\nvxOzommvzIt+eZkVcbMyK3x6ZVb0yzuYWbGN/kP8e2WZF8yKplnDzIpphoXOi1DvxCilfNfM3ifp\nc1p9XscnSilfMbP3rv5y+Xgp5bfM7KfM7GuSntHqXaaQWSVdL+n7Jf0f07uLz5VSLg6adcP/ZdkZ\npZmv/1fN7P+R9HuSvivp46WUP4iYVdIvS7pt3bcT+vullP+07KySZGa/KWlF0n9rZn8o6UatDr1Q\n+2pWzIqmWTf8X5ad8fRvzLxwwaxgViww64b/y7Iznv6NmRUuhjYrJOZFw5wb/i/LzLfhN2ZWuPGY\nFza94wEAAAAAABBatO9OAgAAAAAAsCVuYgAAAAAAgBS4iQEAAAAAAFLgJgYAAAAAAEiBmxgAAAAA\nACAFbmIAAAAAAIAUuIkBF2b222b2Z2Z2b+ssAOJiVgCYBbMCwCyYFePATQx4+SeSDrQOASA8ZgWA\nWTArAMyCWTEC3MTAXMzsIjP7spm9yMzOMrPfN7NXl1KOSHq6dT4AMTArAMyCWQFgFsyKcTujdQDk\nVkr5opndI+kfSXqJpE+WUv6gcSwAwTArAMyCWQFgFsyKceMmBhbhH0p6WNJfSLqmcRYAcTErAMyC\nWQFgFsyKkeI4CRbh5ZJeKmmnpBc3zgIgLmYFgFkwKwDMglkxUtzEwCL8mqR/IOlTWn2Yzik2/QAA\niVkBYDbMCgCzYFaMFMdJMBcze4ek/1pKucPMdkh6wMxWJP2SpB+R9FIz+0NJ/2Mp5b6GUQE0xKwA\nMAtmBYBZMCvGzUoprTMAAAAAAAC8II6TAAAAAACAFLiJAQAAAAAAUuAmBgAAAAAASIGbGAAAAAAA\nIAVuYgAAAAAAgBS4iQEAAAAAAFLgJgYAAAAAAEjh/wdLyO3axv6hXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d6d6a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P1 = mesh._getFacePxx()('fXm','fYm')\n",
    "P2 = mesh._getFacePxx()('fXp','fYm')\n",
    "P3 = mesh._getFacePxx()('fXm','fYp')\n",
    "P4 = mesh._getFacePxx()('fXp','fYp')\n",
    "fig, ax = plt.subplots(1,4, figsize=(15,6))\n",
    "def plot_projection(ii, ax, P):\n",
    "    x = P * np.r_[mesh.gridFx[:,0], mesh.gridFy[:,0]]\n",
    "    y = P * np.r_[mesh.gridFx[:,1], mesh.gridFy[:,1]]\n",
    "    xx, xy = x[:mesh.nC], x[mesh.nC:]\n",
    "    yx, yy = y[:mesh.nC], y[mesh.nC:]\n",
    "    ax.plot(np.c_[xx, mesh.gridCC[:,0], xx*np.nan].flatten(), np.c_[yx, mesh.gridCC[:,1], yx*np.nan].flatten(), 'k-')\n",
    "    ax.plot(np.c_[xy, mesh.gridCC[:,0], xy*np.nan].flatten(), np.c_[yy, mesh.gridCC[:,1], yy*np.nan].flatten(), 'k-')\n",
    "    ax.plot(xx, yx, 'g>')\n",
    "    ax.plot(xy, yy, 'g^')\n",
    "    mesh.plotGrid(ax=ax, centers=True)\n",
    "    ax.set_title('P{}'.format(ii+1))\n",
    "map(plot_projection, range(4), ax.flatten(), (P1, P2, P3, P4));\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each cell (the center is the red dot) has two faces that it uses, one x component and one y component."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The four $P$ matrices are assembled and summed in the **getFaceInnerProduct** code using this formula:\n",
    "\n",
    "\n",
    "$$\n",
    "\\mathbf{M}_f(\\sigma) = \\frac{1}{4}\\sum_{i=1}^{4} \\mathbf{P}_i^\\top \\sqrt{v} \\boldsymbol{\\Sigma}^{-1} \\sqrt{v} \\mathbf{P}_i\n",
    "$$\n",
    "\n",
    "This function takes a sigma variable which can have any sort of anisotropy:\n",
    "\n",
    "$$\n",
    "\\begin{align}\\begin{aligned}\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & 0 \\\\ 0 & \\mu_{1} \\end{matrix}\\right]\\end{split}\\\\\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & 0 \\\\ 0 & \\mu_{2} \\end{matrix}\\right]\\end{split}\\\\\\begin{split}\\vec{\\mu} = \\left[\\begin{matrix} \\mu_{1} & \\mu_{3} \\\\ \\mu_{3} & \\mu_{2} \\end{matrix}\\right]\\end{split}\\end{aligned}\\end{align}\n",
    "$$\n",
    "\n",
    "The matrices that are produced by this function can be seen below!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2QAAADwCAYAAABv96cAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XvYJGV95//3Z2YAORgOrkKYkUFFwRhgwDCLQuRh0XXU\nXXA1KmA8kGjYXUEiufxBTMw8rNlV3GhiQrjCKEEwGERQIW5APD0aNMgAw9EZQJEzjAKKIoLD8P39\nUfXM9PTT1U93V3d13V2f13X1Nd116PuuYfpD3XX4liICMzMzMzMzq96CcXfAzMzMzMysqTwgMzMz\nMzMzGxMPyMzMzMzMzMbEAzIzMzMzM7Mx8YDMzMzMzMxsTDwgMzMzMzMzG5PKBmSSVkhaJ+k2SadU\n1W5L+3dKukHSGklXj7itsyWtl3Rjy7SdJV0h6VZJX5G0Y8Xtr5R0r6Tr8teKEba/RNI3JN0i6SZJ\n782nj/zvoEPbJ+bTK9l+SdtI+l7+7+wmSSvz6ZX997f+NCmb8vYam0/jzKaC9ivLJ2dTepxNzqYm\nZFPejvMpIkb+Ihv4/QBYCmwFXA/sU0XbLX24A9i5orYOBZYBN7ZMOx34//L3pwAfqbj9lcDJFW3/\nbsCy/P0OwK3APlX8HXRpu8rt3y7/cyFwFbC8zLbvCEHvrzur2MZJeTUtm/L2GptP48ymedqvavuH\nmk0RzqcR/rdyNoWzqSnZlLfb6H2nRVRjOXB7RNwFIOkC4ChgXUXtA4iKzghGxJWSlrZNPgo4LH9/\nLjADnFph+5D9HYxcRDwIPJi/f0zSWmAJFfwdFLS9OJ9d1fY/nr/dBlhE9mMfeNsfBf6yx7b/PPuf\nt/WuUdkEzc6ncWZTl/Yry6dhZxM4n0bI2ZRxNjUgm/J2G73vVNUPbTFwT8vne9n8H7oqAXxV0mpJ\n7664bYDnRMR62PQP/zlj6MMJkq6X9KmqTvtK2pPsiNNVwK5V/h20tP29fFIl2y9pgaQ1ZOH21YhY\nTclt36rHl/XN2ZRpXD6NM5va2q8sn0aRTeB8GhFnU8bZ1IBsyttt9L5Tk4p6HBIRBwKvBd4j6dAx\n9ycqbu9M4PkRsYzsH/vHR92gpB2Ai4CTIuIx5m7zyP4OOrRd2fZHxNMRcQDZ0a3lkl5CyW1f1OPL\nklS3bIIJz6dxZlNB+5Vs/yiyCZxPE8zZ5GzyvlNFqhqQ3Qfs0fJ5ST6tMhHxQP7nT4Avkl0OUKX1\nknYFkLQb8OMqG4+In0R+ES7wSeCgUbYnaRHZj/ozEXFJPrmSv4NObVe9/XmbPyc7vb6Cktue0lGe\nxDibMo3Jp3FmU1H7VefTMLMJnE8j4mzKOJsalE15m43cd6pqQLYa2EvSUklbA0cDl1bUNpK2y0f9\nSNoe+M/AzaNuli2vu70UeGf+/h3AJe0rjLL9/B/yrDcw+u3/R+D7EfGJlmlV/R3Mabuq7Zf0H2ZP\n6UvaFngVsJaS275tjy/rWxOzCZqdT+PMpo7tV7H9o8omKJdPmqeSoKQj1VLpT9Ih+fSOVekmiLMp\n42ya8GzK26nlvlOV+aTNA9/RUlYq8xNkg8CzI+IjlTSctf08sqM7QXZ28vxRti/ps8AU8CxgPVmV\nmi8BnweeC9wFvDkiflZh+4eTXRP8NHAncPzsdbkjaP8Q4NvATWyuYPMB4GrgQkb4d9Cl7WOpYPsl\n7Ut24+mC/PW5iPjfknZhwG2XFOf02P5xQERUcgPupGhSNuVtNjafxplN87Q/8nwaRTbl3ztwPkla\nANwGHAHcTzYIOToi1rUss93szf75NlwYES/OdxR3i4jr84HDtcBRreumztnkbKIB2ZS3X7t9p6rz\nqbIBmZkNRlL8U4/L/j4ekJlZdcrkk6SDgZUR8Zr886nZInF6QVsvAz4VES/pMO9LwN9FxNf73QYz\nmzxl952qzqe63MtmZl34h2pmdVUinzpVEpxzn5Kk1wMfBp4NvK7D/D3ZsiqcmVnZfadK86lJVRbN\nkpXSjalm1iyjzqeI+FJEvBh4PW2PFtLcqnBmZkA1+07DyicfeDdLgAdbZlZXRfl0U/7qoq9KgpE9\nuPf5knaJiEcKqtKZmQGlsgkqzqdSZ8jmqz5iZsOR0rM06sL5ZFaNojw6AHh7y6uDeSsJSnpBy/sD\nga0j4pF8UqeqdLXnbDKrRolsgorzaeB9uLz6yBm0VB+RdMkkVTgyqwuXtO+P88msOoPmU0RslHQC\ncAWbKwmulXR8NjtWAW+U9Hbg18CvgDfDpqpwbwVukrSGvCpcRFxecnNGytlkVp0y+05V59PAVRZ7\nrT4iyWUczdr0UwlRUlzZ47KHFnx3Xj75b9gcKu2/0yOBD5GVtt0AvC8iviNpCXAesGs+75MR8be9\n9n1cesknZ5PZXP1WaR1GPjWJ953MBlf1vlOVylyy2Kn6yOLOi67MX4dteh8Rlb9WrhxPu3Vov8nb\nXrf2B1HmksWWI7KvBl4CHCNpn7bFvhYR+0fEAcAfAp/Kpz8FnBxZGdeXAe/psG4d9ZRPh8Vlm15L\nV76Vw+Kysf/7cPvNar9O2z4oX1Ldl573ndqzqYn5VKffR/YbuaHw1e939/Jdddv+cbY/iJSyqS79\nMLMuShb1WA7cHhF3AUi6ADgK2HSJTOQPNsztQHY2jIh4EHgwf/+YpLVkOw++vMbMABcdMrN6Simb\nygzI+qg+MpP/eWf+2rNEs2ZpmZmZYWZmhunp6YG/I6VnadRET/l05/Tmx0Y+cef60ffKrEaGkU3g\nI7t96nnfaTaffjZzIz+buZGdpvYbfe/MaqIG+06VKtPXTdVHgAfIqo8c03nRqfzPOxnnYGxqamre\nZSa1/SZv+7jbn5qaYnp6elMfTjvttL6/o4qjPBHxJeBLkg4le5bGq2bnJfisn57yac/p39/0/mcz\nN1bWuXZN/n00vf1xt102myCto9A10PO+02w+jXswNu5/o+Pk9sfXfir7TsMycFEP2FQo4BNsLhTw\nkQ7LRKc2plV879x0iT6Z1Z0kos8bU28rmPc9tjxddQZzb0zNbyKfjogV+eeON5G3rfND4KDY/CyN\nLwOXRULlpefLJ2eT2Zb6zaZ8ncJ8avcixn/jfB143yldUv8H7iI6D6aH+V1NMMx9p3Z1yKZSZ/Mi\nK9+495D6YmYFtiuYfnj+mnVG58XmPSIr6QUR8cP8/UQ868f5ZFaNonyyzpxNZtVIKZtSurzSrLHK\n/FCjgc/6MbPqeEfCzOoopWxKqa9mjbVVr7/UpzpP7nRENiLOann/UeCjHdb7DrCw136aWfOUzScz\ns1FIKZs8IDNLwKKEQsXMmsX5ZGZ1lFI2eUBmloCtfI7KzGrK+WRmdZRSNnlAZpaAno/ymJlVzPlk\nZnWUUjaVKnvfUwMFpVu7KSrr6pKuNgkGKd0az+px2YfHX7o1Fc4msy0NWvbe+TR8zqf09FvGvlsJ\ne5fEn2vS950SGjuaNdg24+6AmVkB55OZ1VFC2eQBmVkK/Es1s7pyPplZHSWUTQl11azB/Es1s7py\nPplZHSWUTQl11azBEqoUZGYN43wyszpKKJs8IDNLgX+pZlZXziczq6OEsqmWVRaLv2u6cF5E8Tyz\nOhmoUtC+PS570/grBaVimNk0xeWF82ZYMZQ2zEZt4CqLzqehcz5NjqKKiYNURexWfdFVFucsn1Q2\nJTR2NGuwhCoFmVnDOJ/MrI4SyiYPyMxS4F+qmdWV88nM6iihbEqoq2YNltCNqWbWMM4nM6ujhLLJ\nAzKzFPiXamZ15XwyszpKKJsWjLsDZtaDRT2+zMyqViKfJK2QtE7SbZJO6TD/WEk35K8rJe3XMu99\nkm6WdKOk8yVtPfRtM7N0ldx3qjKfPCAzS4EHZGZWVwPmk6QFwBnAq4GXAMdI2qdtsTuAV0TE/sBf\nAqvydXcHTgQOjKy83CLg6CFvmZmlrNzBokrzKalduG6l7YtK4rscvk2EhK6DbqJupaOLSk673LRN\njMHzaTlwe0TcBSDpAuAoYN3sAhFxVcvyVwGL21reXtLTwHbA/QP3ZII5n8arqBz9ICXsu5W2H2Z5\n/YlRbt+p0nxKakBm1ljPGHcHzMwKDJ5Pi4F7Wj7fS7YTVORdwGUAEXG/pI8BdwOPA1dExNcG7omZ\nTZ5y+06V5pMvWTRLwcIeXwV6uA76yPwa6DWSrpZ0SK/rmlnDlcynXkg6HDgOOCX/vBPZ0eqlwO7A\nDpKOLdeKmU2UCrIJhpNPPkNmloISv9SW66CPIDtlvlrSJRGxrmWxr0XEpfny+wIXAi/ucV0za7KC\nfJpZDzM/7rrmfcAeLZ+X5NO2kN8ovwpYERE/zSe/ErgjIh7Jl/kC8HLgs/113swm1uDZBBXnkwdk\nZiko90vt5Trox1uW3wF4utd1zazhCvJpanH2mnXazXMWWQ3sJWkp8ADZTe/HtC4gaQ/gYuBtEfHD\nlll3AwdLegbwJNlBo9WDb4SZTZzBswkqzqdSu3mS7gQeJdt52xAR3a6tNLNBlTul3tN10JJeD3wY\neDbwun7WrSPnk1lFBsyniNgo6QTgCrJbKM6OiLWSjs9mxyrgg8AuwJmSRP5bjoirJV0ErAE25H+u\nKr8xo+dsMqtIiX2nqvNJETFwZyXdAby05RRdp2WiTBtlTUudp4+xT9ZskoiIzv8wOy8f8a7O82bu\nh5kHNn8+bQ1zvlvSG4FXR8Qf5Z9/H1geEe8taO9QYGVEvKrfdetkvnxyNpltqd9sytcpzKc5y35q\nbj41kfedDLpXWSwySJXFQb6rjoa57zRn2RpkU9lLFoULg5iNXtFp9z2y16zT1nRcrKfroGdFxJWS\nni9pl37XrRnnk1kVfPNDv5xNZlVIKJvKBkIAX5W0WtK7h9EhM+tgmx5fnW26Djp/UvzRwKWtC0h6\nQcv7A4Gt85tR5123xpxPZlUol09N5Gwyq0JC2VR27HhIRDwg6dlk4bI2Iq4cRsfMrEWJX2qP10G/\nUdLbgV8DvwLe3G3dUttSHeeTWRUSOgpdE84msyoklE2luhoRD+R//kTSF8lu9p8TKtPT05veT01N\nMTU1VaZZs6TMzMwwMzNT7ktKhkpEXA7s3TbtrJb3HwU+2uu6Kegln5xN1mRDySZIaqenDrzvZDa/\nOuw7VWngoh6StgMWRMRjkrYnO4J+WkRc0bacb0w1azHQjakf6HHZ/zP+G1ProJd8cjaZbWngoh7O\np55538lmuahHfyZ936nM2HFX4IuSIv+e89sDpQ6KwqMobLqtYzY2CR3lqYna55OzySaG86kftc8m\ncD5VYZDBVbdBV9H3DfO7kpNQNg3c1Yj4EbBsiH0xsyIJhUodOJ/MKuR86pmzyaxCCWVTQl01azD/\nUs2srpxPZlZHCWVTQl01a7CalGU1M5vD+WRmdZRQNnlAZpYC/1LNrK6cT2ZWRwllU0JdNWuwhePu\ngJlZAeeTmdVRQtnkAZlZCvxLNbO6cj6ZWR0llE0JdXW4upVnlaY7To/oPN1s5Br7S22ebtk0xeUd\np8+wYlTdMZuf86kxnE/V6LeEfVXflZyEsimhrpo1WEKn3c2sYZxPZlZHCWWTB2RmKfAv1czqyvlk\nZnWUUDYl1FWzBnvGuDtgZlbA+WRmdZRQNnlAZpaChE67m1nDOJ/MrI4SyiYPyMxS4F+qmdWV88nM\n6iihbEqoq9Upqqbo6os2Nv6lGsXVylzdzMbK+WQ4n6pQVDERiqsmFq0zzO+qrYSyKaGumjVYQqfd\nzaxhnE9mVkcJZZMHZGYp8C/VzOrK+WRmdZRQNi0YdwfMrAfb9PgqIGmFpHWSbpN0Sof5x0q6IX9d\nKWm/lnnvk3SzpBslnS9p66Fum5mlrUQ+DZBN+7bNXyDpOkmXDnWbzCx91e87DZxPHpCZpWBRj68O\nJC0AzgBeDbwEOEbSPm2L3QG8IiL2B/4SWJWvuztwInBgZBePLwKOHuKWmVnqBsynAbPpk23zTwK+\nP6QtMbNJUv2+08D55AGZWQpKhAqwHLg9Iu6KiA3ABcBRrQtExFUR8Wj+8SpgccvshcD2khYB2wH3\nD2GLzGxSDJ5PpbJJ0hLgtcCnhrg1ZjYpxrjv1G8+JXR15fgVVVOclgrXmY4YUW+sUcr9UhcD97R8\nvpcsaIq8C7gMICLul/Qx4G7gceCKiPhaqd7Y0BVVK3M2WSUGz6eBsyn318D7gR0H7oGNnPNpvIZZ\nMTG56otj2nfK9ZVPPkNmloKFPb5KknQ4cBxwSv55J7IjQkuB3YEdJB1bviUzmxgV5FOHbHodsD4i\nrgeUv8zMNhvfvlPf+eQzZGYpKPilzlydveZxH7BHy+cl+bQt5IU8VgErIuKn+eRXAndExCP5Ml8A\nXg58tvfOm9lEGzyfymTTIcCRkl4LbAs8U9J5EfH2frtvZhNqfPtOfeeTYsSnhSXFqNsYN592t35I\nIiJ6PporKeL2Hpd9IXO+W9JC4FbgCOAB4GrgmIhY27LMHsDXgbdFxFUt05cDZwMHAU8C5wCrI+Lv\ne+1/XTmbJnvbrX/9ZlO+zsD5VCab2vpwGPAnEXFkP32vM+fTZG/7KBRdTlhkkAdDD/Jdw5LSvlPb\n9/SUTz5DZpaCLmVZ5xMRGyWdAFxBdpny2RGxVtLx2exYBXwQ2AU4U5KADRGxPCKulnQRsAbYkP+5\nqtzGmNlEGTCfymTTcDpuZhNtTPtOg7TnM2RD4KM81o+BjvL0WNdQu889ymOdOZsme9utfwOfIXM+\nDZ3zabK3fRR8hmzO8kll07xFPSSdLWm9Wv7rSNpZ0hWSbpX0FUmucGQ2SuVKt04s55NZDTif5nA2\nmdVAQtk07xkySYcCjwHn5Q+GRdLpwMMR8dH8ydU7R8SpBetP/FGeboqOAPnoT3MNcpTn6Yd7W3bB\ns8Z/lKdKZfLJ2eRssi0NeobM+TSX953KcT4NT79nu6D4jNcwv6tfk77vNO8Zsoi4Evhp2+SjgHPz\n9+cCrx9yv8ysxcZFvb2axvlkNn7Op7mcTWbjl1I2DdqN50TEeoCIeFDSc4bYJzNrU5fASITzyaxC\nzqeeOZvMKpRSNg2rqz6HbDZCTy3s9RnuT4+0H4lyPpmNkPNpYM4msxFKKZsGHZCtl7RrRKyXtBvw\n424LT09Pb3o/NTXF1NTUgM2apWdmZoaZmZlS3/HrbXqt3fqrUu1MiJ7zydlkTTaMbALnUx+872TW\no6btO/VU9l7SnsC/RMS++efTgUci4nTfmNqdb0y1doPcmPqT2KGnZZ+tx8Z+Y2rVBs0nZ5OzybY0\naFEP51Nn3ncanPNpeJpc1COlbOqlyuJngSngWcB6YCXwJeDzwHOBu4A3R8TPCtZvdKgUkaYL50UU\nz7P0DRIqD0Rv1ZF/U4+OPVSqVCafnE2dTXF54bwZVlTYE6vaoAMy59Nc3ncaDefT8HQbXPU7iBrm\ndxW3Mdn7TvNeshgRxxbMeuWQ+2JmBTbW5UEZNeN8Mhs/59Ncziaz8Uspm9LpqVmDbWThuLtgZtaR\n88nM6iilbPKAzCwBKYWKmTWL88nM6iilbPKAzCwBKYWKmTWL88nM6iilbPKAzCwBT9Jr6VYzs2o5\nn8ysjlLKJg/IzBKQ0lEeM2sW55OZ1VFK2eQB2Zh0K21fVBLf5fCbK6VQsbR1Kx1dVHLa5aabzflk\nVXE+DU+3cvRFZeyL1hnmdw1TStnkAZlZAp5KKFTMrFmcT2ZWRyllkwdkZglI6VkaZtYsziczq6OU\nsimdnpo1WEqn3c2sWZxPZlZHKWXTgnF3wMzmt5GFPb2KSFohaZ2k2ySd0mH+sZJuyF9XStqvZd6O\nkj4vaa2kWyT9xxFtppklqGw+mZmNQkrZ5DNkZgl4kq0HXlfSAuAM4AjgfmC1pEsiYl3LYncAr4iI\nRyWtAFYBB+fzPgH8a0S8SdIiYLuBO2NmE6dMPpmZjUpK2eQBWQ0VVVOcljpPjxhhb6wOSl4HvRy4\nPSLuApB0AXAUsGlAFhFXtSx/FbA4X/Y3gN+NiHfmyz0F/LxMZyxdRdXKnE3NltJ9Gja5nE+jN8yK\niUXfNej3dZJSNqXTU7MGK3lKfTFwT8vne8kGaUXeBVyWv38e8JCkc4D9gWuAkyLiV2U6ZGaToy6X\n/JiZtUopm3wPmVkCqroOWtLhwHHA7H1mi4ADgb+PiAOBx4FTSzdkZhOjTD6VvL+167pm1mxjvv++\nr3zyGTKzBBQ9S+PmmYe5ZeaR+Va/D9ij5fOSfNoW8iBZBayIiJ/mk+8F7omIa/LPF7F5sGZmNvCz\nfsrc39rjumbWYGWeQ1Z1PnlAZpaAouugXzy1Ky+e2nXT5wtP+0GnxVYDe0laCjwAHA0c07qApD2A\ni4G3RcQPZ6dHxHpJ90h6UUTcRhYu3y+3NWY2SUrcpzHw/a29rGtmzTau++97WbedB2RmCfh1iUpB\nEbFR0gnAFWSXKZ8dEWslHZ/NjlXAB4FdgDMlCdgQEbP3mb0XOF/SVmRHg44rsSlmNmFK5FOZ+1v7\nXdfMGqbMvhMV55MHZAkpqghUVEGo2zqWljKn3QEi4nJg77ZpZ7W8fzfw7oJ1bwAOKtUBm2jOpmYr\nm0+9aLm/9dCRN2YTxfnUv6Iqh0WVEQepmNhtnW7z+lFFNsFw8skDMrMEpFS61cyapSif1s2s59aZ\n9d1WLXN/a0/rmllzlcgmqDifvJdnloCUSreaWbMU5dMLp3bnhVO7b/p86Wk3ty8y8P2tvaxrZs1W\nIpug4nzygMwsAR6QmVldDZpPZe5vLVp3GNtjZpOhzL5T1fnkAZlZAjwgM7O6KrnTU+b+1jnrmpnN\nKrvvVGU+eUBmloCqbkw1M+uX88nM6iilbFow3wKSzpa0Xi0lTyStlHSvpOvy14rRdtOs2X7NNj29\nmsb5ZDZ+zqe5nE1m45dSNvVyhuwc4O+A89qmfzwiPj78Llm/upVnlaY7To/oPN3qyZcsFnI+1Vi3\nbJri8o7TZ/A+amqcTx05m2rO+dS/QUrY9/tdg35fJyll07wDsoi4Mq8S0q74AQ5mNlQpnXavkvPJ\nbPycT3M5m8zGL6VsmveSxS5OkHS9pE9J2nFoPTKzOTayqKeXbeJ8MquI86kvziaziqSUTYP24kzg\nf0VESPpL4OPAHxYtPD09ven91NQUU1NTAzZrlp6ZmRlmZmZKfUdKp91roOd8cjZZkw0jm8D51Afv\nO5n1bDVwTalvSCmbFF2uod20UHba/V+iwwWf3ebl86OXNmw0fA9Z/UgiInq+bEVSfCA+2NOy/0cf\n6uu7J8Gg+eRsGi/fo1E//WZTvo7zqYD3ndLlfOpPt3u+ut0r1t/37T/R+069niETLdc9S9otIh7M\nP74B6PiIazMbjpSO8oyB88lsjJxPhZxNZmOUUjbNOyCT9FlgCniWpLuBlcDhkpYBTwN3AsePsI9W\nQtGZsKIzZ93WsfF5siZlWevG+ZSuoiPNzqb0OJ/mGkY29VtpbpAzEdaZ86k/w/631+n7NMD5q5Sy\nqZcqi8d2mHzOCPpiZgVSOspTJeeT2fg5n+ZyNpmNX0rZVI/SImbWVUqhYmbN4nwyszpKKZs8IDNL\nQErP0jCzZnE+mVkdpZRNHpCZJaAuz8kwM2vnfDKzOkopm9LpqVmDpXTa3cyaxflkZnWUUjYtGHcH\nzGx+G1nY06uIpBWS1km6TdIpHeYfK+mG/HWlpH3b5i+QdJ2kS0eweWaWsLL5ZGY2Cillk8+QNVS3\n8qx+mHT9PMnWA68raQFwBnAEcD+wWtIlEbGuZbE7gFdExKOSVgCfBA5umX8S8H3gNwbuiFkPnE3p\nKZNPVqxz6e/iUvhF81wOf3icT2lJKZs8IDNLQMnroJcDt0fEXQCSLgCOAjYNyCLiqpblrwIWz36Q\ntAR4LfC/gZPLdMTMJk9K92mYWXOklE3p9NSswUqeUl8M3NPy+V6yQVqRdwGXtXz+a+D9wI5lOmFm\nk6kul/yYmbVKKZs8IDNLQFWhIulw4Djg0Pzz64D1EXG9pClAlXTEzJKR0k6PmTVHStnkAZlZAoqe\npfHwzM08MnPLfKvfB+zR8nlJPm0LkvYDVgErIuKn+eRDgCMlvRbYFnimpPMi4u39bYGZTaqUnvVj\nZs2RUjZ5QGaWgKLroHeaWsZOU8s2ff7BaRd2Wmw1sJekpcADwNHAMa0LSNoDuBh4W0T8cHZ6RHwA\n+EC+zGHAn3gwZmatUrpPw8yaI6VsSqenVpmiikCuIDQ+ZU67R8RGSScAV5A96uLsiFgr6fhsdqwC\nPgjsApwpScCGiOh2n5lZ5ZxN9ZTSZUGp61YxsajKYrfKjK7AODzOp/pJKZv8HDKzBDzJ1j29ikTE\n5RGxd0S8MCI+kk87Kx+MERHvjohnRcSBEXFAp8FYRHwrIo4c2UaaWZLK5FMPz0jcW9J3JT0h6eS2\neTtK+ryktZJukfQfR7SJZpagsvtOVeaTz5CZJSCl0+5m1iyD5lOPz0h8GDgReH2Hr/gE8K8R8SZJ\ni4DtBuqImU2kMvtOVeeT9/LMEpDSaXcza5YS+dTLMxIfAh6S9F9aV5T0G8DvRsQ78+WeAn4+aEfM\nbPKU3HeqNJ88IDNLgAdkZlZXJfKp32cktnoe2Y7QOcD+wDXASRHxq0E7Y2aTpeJnuLbqO588IDNL\ngAdkZlZXRfn0i5nreGzmulE1uwg4EHhPRFwj6W+AU4GVo2rQzNIypmyCAfLJAzLrWVFFoGkVPyt4\nOmJEvWmWlJ6lYVa1wupmO3dZ56fF86w/Rfm07dRBbDt10KbP6087u32Rnp6RWOBe4J6IuCb/fBEw\n56b7lHWrjtjJIBUTq2ij6ZxP41Mim6DifPKAzCwBv2abcXfBzKyjEvk07zMS22w6+hcR6yXdI+lF\nEXEb2Y333x+0I2Y2eUruO1WaTx6QmSXAlyyaWV0Nmk+9PCNR0q5k9188E3ha0knAb0XEY8B7gfMl\nbQXcARw3hM0xswkx6me4DjOfPCAzS4AvWTSzuiqTTxFxObB327SzWt6vB55bsO4NwEGd5pmZld13\nqjKfPCAzS4CfQ2ZmdeV8MrM6SimbFsy3gKQlkr6RP2X6JknvzafvLOkKSbdK+oqkHUffXbNm2sjC\nnl5N4mwOq33kAAAaMklEQVQyqwfn01zOJ7PxSymb5h2QAU8BJ0fES4CXAe+RtA9Z+cavRcTewDeA\nPx1dN82aLaVQqZCzyawGnE8dOZ/MxiylbFL0WZZc0peAM/LXYXklkd2AmYjYp8Py0W8bNjmKSuI3\nuRy+JCKi+FkBc5ePZ228t6dlH164pK/vniTOJutHUcnpJpeb7jeb8nWcTz0YVj71W6YeqimHP2g7\n1pnzaa5J33fq6+JKSXsCy4CrgF3zm9mIiAclPWfovTMzAJ58wmXvu3E2mY2P86k755PZeKSUTT0P\nyCTtQPZgs5Mi4jFJ7YdufKjZbEQ2PlWPU+p15GwyGy/nUzHnk9n4pJRNPQ3IJC0iC5TPRMQl+eT1\nknZtOe3+46L1p6enN72fmppiampq4A6bpWZmZoaZmZlS35FSqFTJ2WQ2uGFkEzifijifzAbXtH2n\nnu4hk3Qe8FBEnNwy7XTgkYg4XdIpwM4RcWqHdX2fRoP5HrK5BrkOesGDj/W07NO77TD266Cr5Gyy\nQfkejbkGvYfM+dTZKPLJ95A1g/Nprknfd5r3DJmkQ4C3AjdJWkN2ev0DwOnAhZL+ALgLePMoO2rW\nZE9vTOdZGlVxNpnVg/NpLueT2fillE19V1nsuwEfhbYOis6cweSfPRvkKA93beht4aVbjf0oTyqc\nTdZJ0ZFpmPyj04OeIXM+Dd8g+VTFWS2fORsv59Pk7jv18hwyMxu3pxb29iogaYWkdZJuyy+TaZ9/\nrKQb8teVkvbNp3d8uKmZ2SYl88nMbCQSyqZ0zuWZNdkTgx+4kbSA7Nk3RwD3A6slXRIR61oWuwN4\nRUQ8KmkF8EngYDY/3PT6vFrYtZKuaFvXzJqsRD6ZmY1MQtnkAZlZCp4qtfZy4PaIuAtA0gXAUcCm\nQVVEXNWy/FXA4nz6g8CD+fvHJK3N53lAZmaZcvlkZjYaCWWTB2RmKSgXKouBe1o+30s2SCvyLuCy\n9oktDzf9XqnemNlkSWinx8waJKFs8oDMLAU93pdalqTDgeOAQ9umb/Fw02p6Y2ZJqCifzMz6klA2\neUBmloKNBdOvm4E1M/OtfR+wR8vnJfm0LUjaD1gFrIjYXK+p4OGmZmaZonwyMxunhLLJZe+tdib9\nYdIDlW79Vo/bftjc75a0ELiVrKjHA8DVwDERsbZlmT2ArwNva7ufrOPDTSeBs8n6NekPax247H2J\nfLLOqsqnfsvY+yHT9eV8mrN8UtnkM2RmKShxHXREbJR0AnAF2aMuzo6ItZKOz2bHKuCDwC7AmZIE\nbIiI5UUPN42Iy8ttkJlNjITu0zCzBkkomzwgM0vBE+VWzwdQe7dNO6vl/buBd3dY7ztAPR7SYWb1\nVDKfzMxGIqFs8oDMLAUJHeUxs4ZxPplZHSWUTR6QmaUgoVAxs4ZxPplZHSWUTQvG3QEz68GGHl9m\nZlUrkU+SVkhaJ+k2Sad0mL+3pO9KekLSyS3Tl0j6hqRbJN0k6b1D3y4zS1vJfacq88lnyKx2iqop\nTnr1xa4SKt1qNqmKqpVNenWzeQ2YT5IWAGeQVYC9H1gt6ZKIWNey2MPAicDr21Z/Cjg5Iq7Pn5N4\nraQr2ta1eRRVMyyqjNitYmLRd/XbRrd5rr5YzPnUQYl9p6rzyWfIzFLwVI8vM7OqDZ5Py4HbI+Ku\niNgAXAAc1bpARDwUEde2f0NEPBgR1+fvHwPWAouHt1Fmlrxy+06V5pPPkJmlwIMtM6urwfNpMXBP\ny+d7yXaC+iJpT2AZ8L2Be2Jmk6fcvlOl+eQBmVkKEirdamYNM8Z8yi8Hugg4KT8SbWaWGfO+Uz/5\n5AGZWQp8hszM6qoon26fgR/MdFvzPmCPls9L8mk9kbSIbGfnMxFxSa/rmVlDDJ5NUHE+eUBmlgIP\nyMysrory6XlT2WvW5ae1L7Ea2EvSUuAB4GjgmC4ttVd2+kfg+xHxid47a2aNMXg2QcX55AGZJaPf\n6ovd1kmOS9qb1Va/1c26rZOkAfMpIjZKOgG4gqzI2NkRsVbS8dnsWCVpV+Aa4JnA05JOAn4L2B94\nK3CTpDVAAB+IiMtLb49VUhmxW8XEYVZ5bLpG51OJfaeq88kDMrMUuOy9mdVViXzKd1D2bpt2Vsv7\n9cBzO6z6HWDh4C2b2cQrue9UZT55QGaWAl+yaGZ15XwyszpKKJs8IDNLgassmlldOZ/MrI4SyiYP\nyMxS4HvIzKyunE9mVkcJZdOC+RaQtETSNyTdIukmSSfm01dKulfSdflrxei7a9ZQG3t8NYizyawm\nnE9zOJ/MaiChbOrlDNlTwMkRcX3+gLNrJX01n/fxiPj46LpnZkBS10FXyNlkVgfOp06cT2bjllA2\nzTsgi4gHgQfz949JWgsszmcX1xs3q0i30vZFJfGTK4efUKhUxdlkddetdHRRyekky007n+aY1Hwa\npLR8t1L142yn6WXyG5FPCWXTvJcstpK0J7AM+F4+6QRJ10v6lKQdh9w3M5u1ocdXQzmbzMbI+dSV\n88lsTBLKpp4HZPkp94uAkyLiMeBM4PkRsYzsKJBPv5uNSsnroCWtkLRO0m2STukw/1hJN+SvKyXt\n1+u64+ZsMhuzhO7TqJrzyWyMEsqmnqosSlpEFiifiYhLACLiJy2LfBL4l6L1p6enN72fmppiampq\ngK6apWlmZoaZmZlyX1KidKukBcAZwBHA/cBqSZdExLqWxe4AXhERj+Y3ma8CDu5x3bFxNpkNbijZ\nBEmVlq6S88lscOPed6qaood7aSSdBzwUESe3TNstv0YaSe8DDoqIYzusG720YTYKdbyHTBIR0fM9\nBJKC3+uxvxfN/W5JBwMrI+I1+edTgYiI0wva2wm4KSKe2++6VXM2WarqeI9Gv9mUr1MqnyaZ8ynj\ne8jSMwn5lFo2zXuGTNIhwFuBmyStAQL4AHCspGXA08CdwPEj7KdZs5U7pb4YuKfl873A8i7Lvwu4\nbMB1K+NsMquJmlzyUyfOJ7MaSCibeqmy+B1gYYdZlw+/O2bDVXQmrI5nzroqqhT00Aw8PDO0ZiQd\nDhwHHDq0Lx0RZ5OlrOhIcx2PTM8roUpmVXE+bVZ0JqrbGa2ied3OavXbTrf2m372bGLyKaFs6uke\nMjMbs6JQ2Wkqe8267bROS90H7NHyeUk+bQt5IY9VwIqITfHa07pm1mAJ7fSYWYMklE0ekJmloFxZ\n1tXAXpKWAg8ARwPHtC4gaQ/gYuBtEfHDftY1s4arSdloM7MtJJRNHpCZpaDEddARsVHSCcAVZI+6\nODsi1ko6Ppsdq4APArsAZ0oSsCEilhetW3JrzGySJHSfhpk1SELZ5AGZWQpKlm6NiMuBvdumndXy\n/t3Au3td18xsk4RKS5tZgySUTR6QmaUgodPuZtYwziczq6OEsskDMmukfqsvdlunEgmddjezwfVb\n3azbOpVxPtkAulUyHGZlxKqqPDZBcvmUUDZ5QGaWgoQqBZlZwzifzKyOEsomD8jMUpBQqJhZwzif\nzKyOEsomD8jMUpDQddBm1jDOJzOro4SyyQMysxQkdB20mTWM88nM6iihbFow7g6YWQ+e6PFlZla1\nEvkkaYWkdZJuk3RKwTJ/K+l2SddLWtYy/X2SbpZ0o6TzJW091O0ys7SV3HeqMp88IDNLwYYeX2Zm\nVRswnyQtAM4AXg28BDhG0j5ty7wGeEFEvBA4HviHfPruwInAgZGVvlsEHD38jTOzZJXYd6o6n3zJ\nolmLbqXtpemO0yM6Tx+qhE67m9nwdSsdPdZsgjL5tBy4PSLuApB0AXAUsK5lmaOA8wAi4nuSdpS0\naz5vIbC9pKeB7YD7B+6J1cog5eW7lbEfVhvWWW3zqdy+U6X55DNkZil4qseXmVnVBs+nxcA9LZ/v\nzad1W+Y+YHFE3A98DLg7n/aziPhayS0xs0lSbt+p0nzygMwsBR6QmVldjSGfJO1EdnR6KbA7sIOk\nY4fbipklbUz7ToPkky9ZNEuB7w8zs7oqyqenZyBmuq15H7BHy+cl+bT2ZZ7bYZlXAndExCMAkr4A\nvBz4bK/dNrMJN3g2QcX55AGZWQp8D5mZ1VVhPk3lr1mntS+wGthL0lLgAbKb3o9pW+ZS4D3A5yQd\nTHbpz3pJdwMHS3oG8CRwRP59ZmaZwbMJKs4nD8jMUlBca8TMbLwGzKeI2CjpBOAKslsozo6ItZKO\nz2bHqoj4V0mvlfQD4JfAcfm6V0u6CFhDdhx8DbCq/MaY2cQose9UdT4pulSVGwZJMeo2zMapqIIQ\ndK4iJImIUO/fr+g9Vfr77iZzNtmkG3U25es4n0bA+TQ8/VZfBFdgrIL3nbbkoh5mZmZmZmZj4ksW\nzZLgqh5mVlfOJzOro3SyyWfIzJJQrnarpBWS1km6TdIpHebvLem7kp6QdHLbvPdJulnSjZLOl7T1\n8LbLzNLn53KYWR2lk00+Q2aWhMGP8khaAJxBVuXnfmC1pEsiovVp8w8DJwKvb1t393z6PhHxa0mf\nI6s0dN7AHTKzCZPOUWgza5J0smneM2SStpH0PUlrJN0kaWU+fWdJV0i6VdJXJO04+u6aNVWpozzL\ngdsj4q6I2ABcQPbAwk0i4qGIuLbgSxYC20taBGxHNqgbO2eTWV2kcxS6Ks4nszpIJ5vmHZBFxJPA\n4RFxALAMeI2k5cCpwNciYm/gG8CfjrSnZo32eI+vjhYD97R8vjefNq+IuB/4GHA32cMOfxYRX+u/\n/8PnbDKri1L5NJGcT2Z1kE429XTJYkTM9nabfJ0gO8J+WD79XGCGLGjMGqVTedZZ3cq69mc8R3Ak\n7UT2W18KPApcJOnYiCh82nyVnE1mxarJJqjLEea6cT7VQ1EJ+27l8IvmuRz+8EzyvtMgehqQ5feg\nXAu8APj7iFgtadeIWA8QEQ9Kes4I+2nWcEXXQX8vf3V1H7BHy+cl+bRevBK4IyIeAZD0BeDlQC0G\nZM4mszpI5z6NKjmfzMYtnWzqqcpiRDydn3ZfAiyX9BLmPm1t3qevzczM9N3BYWpy+03e9jq0D3eW\nXL/ouueXAv+z5dXRamAvSUvzColHA5d2aaz14Yh3AwdLeoYkkRUGWVtiQ4ZqWNkE/n26/ea1nblz\nCN+Rzn0aVfK+U9ptZ1aPtfVxb/+42x/dvlP9sqmvsvcR8XOy0+srgPWSdgWQtBvw46L1pqenN72a\n/MP2tje1/TvJfjazr0Fs6PE1V0RsBE4ArgBuAS6IiLWSjpf0RwCSdpV0D/A+4M8k3S1ph4i4GrgI\nWAPcQDZYWzXgRozMMLLp05/+dBVd7ajZv49mt59+NkGZfGoC7zul2XbmmrG2Pu7tTz+f0smmeS9Z\nlPQfgA0R8aikbYFXAR8hO8L+TuB04B3AJUXfMT09venPqampsn02S8ye+Wsq//ytAb6j3BGciLgc\n2Ltt2lkt79cDzy1Y9zTgtFIdGIFhZlP7e7Nm2JPy2QR1OcJcJ953MitrT8a971SlXu4h+03g3Pxa\n6AXA5yLiXyVdBVwo6Q+Au4A3j7CfZg33q3F3oI6cTWa14HzqwPlkNnbpZJMierq9YvAGpNE2YJag\niND8S2Wy39A3e1z68L6+u8mcTWZz9ZsfzqfRcD6ZzTXJ+049VVksY9wbaDYZ0jntngpnk9mwOJ+G\nzflkNgzpZNPIB2RmNgz1uOnUzGwu55OZ1VE62eQBmVkS0jnKY2ZN43wyszpKJ5s8IDNLQjpHecys\naZxPZlZH6WRTX88hs+pI+sWA6/3pkPtxmqT/NMzvtEGk8ywNG438WXH/LOl2SaslfVnSXiNoZ6Wk\nk/P3A//+Je0v6TUDrvvHkn4l6Zk9Lv9lSb8xQDsnSXpG/z20LTmfbDgkbZR0naQ1+Z97zLP8jyTt\nkr8v3G+S9HpJT0t6UY/9WCVpn/56D5LekT9fzmohnWzyGbL6GrTC0geAD3eaIUnRZ1nNiFg5YD9s\nqNIp3Woj80XgnIg4BkDSvsCuwA8G/UJJC/MHh3dU8ve/DPgd4LIB1j0auBp4A3DufAtHxH8ZoA2A\nPwY+AzzRPkPSgoh4esDvbRjnkw3NLyPiwD6Wj4L37Y4G/g04hh6eqxkRf9RHH1q9E7gZeLB9hjNl\nHNLJJp8hqzlJu0n6Vn6k6EZJh+TTj8k/3yjpw/m0DwPb5st+RtJSSesknSvpJmBJ23ofaWnnF5I+\nLulmSV+V9Kx8+jmS3pC/P0jSdyRdL+kqSdtX/zfSVE/1+LJJJOlw4NcR8cnZaRFxU0R8J5//fyXd\nJOkGSW9uWW/OdEmHSfq2pEuAW/JpfybpVknfpuUB4m2//x9JmpZ0bf59L8qnHyTpu/n0KyW9UNJW\nwP8C3pzn0ZskbSfp7Dw7rpX0Xwu29fnA9sCfA8e2TH+HpIslXZb39fSWeT+StEvexpfzo+s3SnpT\nPv+IvB83SPqUpK0lnQjsDnxT0tfz5X4h6a8krQEO7rDeVi3tnZ63cZWk50vaQdIdkhbmyzyz9fNk\ncz7Z0MypLpn/9v+u5fO/SHpF0fId1t8eOAT4Q7IB2ez0wyR9U9LnJa2V9JmWed+UdKCkBXkO3pjn\nwEn5/GWS/j3fH7pY0k6S3kh2EOqf8tx4Rp4VH5F0DfB7yq4caF1vx5b2/qYlu35Hmdta9sek7AqJ\nZw3w99pQ6WSTB2T1dyxweX7EaH/gekm/CXyE7PHly4Dlko6MiD8FHo+IAyPibfn6ewFnRMS+ZP/q\nWtc7SNKR+XLbA1dHxG8D3wa2ODKe74hcAJwYEcuAV5LSoYfkpXPa3Ubit4FrO83IB0z75b/xVwH/\nV9nljR2n56sdQPZb3kfSgWQPp90PeB1wUJd+/DgiXgr8A/D+fNpa4NB8+krgwxGxAfgLsofhHhgR\nnwf+DPh6RBwM/CfgryRt26GNo4F/Bq4EXiTp2S3z9gfelPf1LZIW59Nnj4yvAO6LiAMiYj/gcknb\nAOcAb4qI/YGtgP8eEX8H3A9MRcQR+frbA/8eEQeQ/X23r/c/Wvry07yNvwc+ERGPkT305nUt23Fx\ntzOQk8P5ZEMze1B5jaSLW6aXeS7bUWT7UT8AHpJ0QMu8ZcB7gd8CXiDp5W3rLgMWR8R+eQ6ck08/\nF3h/vj90M/AXEXExcA1wbJ57s2feH4qI34mIC4Hz2tZr3dfaNs+e95BdDRFkZ/B/P5//SuD6iHi4\nxN9Fw6STTR6Q1d9q4DhJf0G2c/VLsh2mb0bEI/np7/OBoqNFd0XE6vx9t/WeBi7M3/8TcGjb9+wN\n3B8R1wFExGM+9V6ldI7yWOUOJRvAEBE/BmaA5QXTZwdbV0fE3fn73wW+GBFPRsQvgEu7tPXF/M9r\ngaX5+52Ai5Sdhf9rsh2bTv4zcGp+9mkG2BrodH/IMWQDuQC+QDYAm/X1PHueBL7f0ofZ3LsJeJWk\nD0s6NN+evYE7IuKH+TLnsjn3WteF7Ef0hfz9fOtdkP/5z8DL8vdnA8fl749j887bhHM+2dDMHlQ+\nICLeOKTvPIbNv9fP0XLmnSwLH8jz5npgz7Z17wCeJ+kTkl4N/ELZ/ao7RsSV+TLdMmW2TXpYbzav\n/w14Zr78OcDsAfY/oDGZMizpZJPvIau5iPi3/NT864BzJH0c+Dk9nKbP/bLtc6/rdToa5QdVjk09\njuDY2NwC/F6Py4r5f7/tudCrJ/M/N7L5/x8fAr4REW+QtJTsLFGRN0bE7UUzJf028ELgq5IgG7T9\nCDizrf32PgAQEbfnZ/xeC3wovxTxUnrPrifa7rPttl7rck/n7X9X0p6SDgMWRMT3e2w3cc4nG6mn\n2PIEQs+FeCTtTHZG/rclBbCQ7Lc7e4Z/vkz5maT9gVcD/53sANHJ9Lc/1GvetmdPRMS9ktYru2z9\nILYcTNq80skmnyGrLwEoqzD044g4m+zo64FkN7u/Ir9nYiHZ0Z+ZfL1ft92z0Boa3dZbwOYdvreS\nXS7U6lZgN0kvzfu1gyT/+6lMOkd5bPgi4hvA1pLeNTtN0r6SDiW7Uf0t+b0OzyY743V1l+ntvg28\nXtI2yqoadry3q4sdgfvy98e1TP8F0Fr58CtklwbN9n9Zh+86BlgZEc/PX0uA3SU9t5eO5Jdz/yoi\nPgv8FVle3gosVXZvGmRHm2fy9z9v62NrXnZbD+At+Z9HA//eMv0zwGeBf+ylz5PB+WRD02mgcyew\nLL+H6rlkVwD0uu6bgPMi4nl5piwFfpRn5/ydye7XWhgRXyS7r/XAiPg58Ijye/rJsuFb+fv23Nsk\nX++nBetBnil5336Wn+GHbN/vn4AL+y3MZulkk8+Q1dfsj24KeL+kDWQ/9LdHxIOSTmXzzsGXI+LL\n+ftVwE2SriULj00/3g7r/b+W9X5Jdi/aB4H1bN7ZiHzdDZLeApyR3/fxONn1zI8Pb5OtWDpHeWxk\n/hvwifw3/CuynZQ/jogrJb0MuIHsTM3780sUvyjp4Pbpkl7c+qURsUbS54AbyX77rYO2XiqYfRQ4\nV9KfA/+vZfo3yS5RvI6s8uuH8v7fSLbj9CPgyLbvegvZ2a1WXyQb9Kxvm96pb/uS3Sv3NPBr4H9E\nxJOSjiO7rHIh2WXgZ+XLf5LsPrP78vvIWvOy23oAO0u6gaxC4zEt08/Pt/UCGsP5ZEMzJ2ci4juS\n7iS7UmAtW95PO19GvQU4vW3axWS/2Qvbpnf6rsVkVyctyKedmk9/J/AP+f7QHWw+GPXpfPrjwMs7\n9OkdwFkd1gN4Is/LRW3TLyU7wPPpDttnXaWTTfJg2yCrLhYRPT3zx6qV/49o6XzL5e6KiD1H1xsz\nk/Qj4KUR8UiHeb8H/NeIeEf1Paue88msPEnfBP5k9j79tnm/A3wsIg6rvmfpSi2bfIbMZnlkXlPj\nDgkzm6NjXkr6W7JKj+1n+SaW88lsKIoy5RSye9d871ifUssmnyEzMzMzMzMbExdlMDMzMzMzGxMP\nyMzMzMzMzMbEAzIzMzMzM7Mx8YDMzMzMzMxsTDwgMzMzMzMzGxMPyMzMzMzMzMbk/wc6CBDX7Atn\nIQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fa5b950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "isotropic       = np.ones(mesh.nC)*4\n",
    "anisotropic_vec = np.r_[np.ones(mesh.nC)*4, np.ones(mesh.nC)]\n",
    "anisotropic     = np.r_[np.ones(mesh.nC)*4, np.ones(mesh.nC), np.ones(mesh.nC)*3]\n",
    "Mf_sig_i = mesh.getFaceInnerProduct(isotropic)\n",
    "Mf_sig_v = mesh.getFaceInnerProduct(anisotropic_vec)\n",
    "Mf_sig_a = mesh.getFaceInnerProduct(anisotropic)\n",
    "clim = (\n",
    "    np.min(map(np.min, (Mf_sig_i.data, Mf_sig_v.data, Mf_sig_a.data))),\n",
    "    np.max(map(np.max, (Mf_sig_i.data, Mf_sig_v.data, Mf_sig_a.data)))\n",
    ")\n",
    "fig, ax = plt.subplots(1,3, figsize=(15,4))\n",
    "def plot_spy(ax, Mf, title):\n",
    "    dense = Mf.toarray()\n",
    "    dense[dense == 0] = np.nan\n",
    "    ms = ax.matshow(dense)#, clim=clim)\n",
    "    plt.colorbar(ms, ax=ax)\n",
    "    ax.set_xlabel(title)\n",
    "map(plot_spy, ax, (Mf_sig_i, Mf_sig_v, Mf_sig_a), ('Isotropic', 'Coordinate Anisotropy', 'Full Anisotropy'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last figure you can see that there is coupling between the cells in x and y due to the non-coordinate aligned anisotropy. This will make the matrix system harder to solve when we try to invert it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ok, but does it work?\n",
    "\n",
    "Unfortunately the testing of the innerproducts is not that visual. We need to make sure that the integral of some made up function evaluates to the correct value.\n",
    "\n",
    "So to make something up, and integrate it, we will use **sympy**!\n",
    "\n",
    "We will evaluate the integral:\n",
    "$$\n",
    "\\int_0^1 \\int_0^1 \\sigma(x,y) \\vec{j}(x,y) \\cdot \\vec{j}(x,y) \\quad dx dy\n",
    "$$\n",
    "by making up some ridiculous functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It is trivial to see that the answer is 42.\n"
     ]
    }
   ],
   "source": [
    "import sympy\n",
    "from sympy.abc import x, y\n",
    "\n",
    "# Here we will make up some j vectors that vary in space\n",
    "j = sympy.Matrix([\n",
    "    x**2+y*5,\n",
    "    (5**2)*x+y*5\n",
    "])\n",
    "\n",
    "# Create an isotropic sigma vector\n",
    "Sig = sympy.Matrix([\n",
    "    [x*y*432/1163,      0      ],\n",
    "    [     0      , x*y*432/1163]\n",
    "])\n",
    "\n",
    "# Do the inner product!\n",
    "jTSj = j.T*Sig*j\n",
    "ans  = sympy.integrate(sympy.integrate(jTSj, (x,0,1)), (y,0,1))[0] # The `[0]` is to make it an int.\n",
    "\n",
    "print( \"It is trivial to see that the answer is {}.\".format(ans) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wonderful. Everything seems right in the world with that answer. Next step is to sub in the locations of the grid on both the faces in x and y as well as the cell centers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_vectors(mesh):\n",
    "    \"\"\"Gets the vectors sig and [jx, jy] from sympy.\"\"\"\n",
    "    f_jx  = sympy.lambdify((x,y), j[0], 'numpy')\n",
    "    f_jy  = sympy.lambdify((x,y), j[1], 'numpy')\n",
    "    f_sig = sympy.lambdify((x,y), Sig[0], 'numpy')\n",
    "    jx  = f_jx(mesh.gridFx[:,0], mesh.gridFx[:,1])\n",
    "    jy  = f_jy(mesh.gridFy[:,0], mesh.gridFy[:,1])\n",
    "    sig = f_sig(mesh.gridCC[:,0], mesh.gridCC[:,1])\n",
    "    return sig, np.r_[jx, jy]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these vectors ($\\sigma$ and $\\mathbf{j}$) we can evaluate the result on our mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerically we get 41.1891755804.\n"
     ]
    }
   ],
   "source": [
    "n = 5 # get's better if you add cells!\n",
    "mesh = Mesh.TensorMesh([n,n])\n",
    "sig, jv = get_vectors(mesh)\n",
    "Msig = mesh.getFaceInnerProduct(sig)\n",
    "numeric_ans = jv.T.dot(Msig.dot(jv))\n",
    "print( \"Numerically we get {}.\".format(numeric_ans) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is pretty close, but to be rigorous we should look at the order of the convergence which should be $\\mathcal{O}(h^2)$.\n",
    "\n",
    "SimPEG has a number of testing functions for \n",
    "[derivatives](http://docs.simpeg.xyz/content/api_core/api_Tests.html#SimPEG.Tests.checkDerivative)\n",
    "and \n",
    "[order of convergence](http://docs.simpeg.xyz/content/api_core/api_Tests.html#SimPEG.Tests.OrderTest) \n",
    "that make this pretty simple!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "."
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "uniformTensorMesh:  Order Test\n",
      "_____________________________________________\n",
      "   h  |    error    | e(i-1)/e(i) |  order\n",
      "~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~\n",
      "   4  |  1.27e+00   |\n",
      "   8  |  3.17e-01   |   3.9942    |  1.9979\n",
      "  16  |  7.93e-02   |   3.9985    |  1.9995\n",
      "  32  |  1.98e-02   |   3.9996    |  1.9999\n",
      "---------------------------------------------\n",
      "Awesome, Rowan, just awesome.\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "----------------------------------------------------------------------\n",
      "Ran 1 test in 0.028s\n",
      "\n",
      "OK\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import unittest\n",
    "from SimPEG.Tests import OrderTest\n",
    "\n",
    "class Testify(OrderTest):\n",
    "    meshDimension = 2\n",
    "    def getError(self):\n",
    "        sig, jv = get_vectors(self.M)\n",
    "        Msig = self.M.getFaceInnerProduct(sig)\n",
    "        return float(ans) - jv.T.dot(Msig.dot(jv))\n",
    "    def test_order(self):\n",
    "        self.orderTest()\n",
    "\n",
    "# This just runs the unittest:\n",
    "suite = unittest.TestLoader().loadTestsFromTestCase( Testify )\n",
    "unittest.TextTestRunner().run( suite );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These basic operations on the mesh for inner products and differential operators are core pieces of SimPEG, and they need to be tested continuously to ensure that no one has broken them!\n",
    "\n",
    "On every single change to the SimPEG codebase many hours of testing are completed!\n",
    "\n",
    "If you would like to browse through the results of 424 tests that are run on the mesh classes you can see them here:\n",
    "\n",
    "https://travis-ci.org/simpeg/simpeg/jobs/148608973"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next up ...\n",
    "\n",
    "In the [next notebook](all_together_now.ipynb) we will use our knowledge of the [mesh](mesh.ipynb) and the [divergence](divergence.ipynb) implementations to bring the DC resistivity equations together and show some pretty figures!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
